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1.  INTRODUCTION  AND  RESEARCH  OBJECTIVES 


This  report  describes  work  carried  out  at  the  Gas  Turbine  and  Plasma 
Dynamics  Laboratory  at  MIT,  as  part  of  our  multi-investigator  effort  on 
current  problems  in  turbomachinery  fluid  dynamics.  Support  for  this 
program  is  provided  by  the  Air  Force  Office  of  Scientific  Research  under 
Contract  Number  F49620-82-K-0002,  Major  M.S.  Francis  and  Dr.  J.  D.  Wilson, 
Program  Managers. 

The  present  report  gives  a  summary  of  the  work  for  the  period  11/1/83 
4/30/84.  For  further  details  and  background,  the  referenced  reports,  pub¬ 
lications  and  previous  reports  (covering  the  period  10/1/79  -  5/31/83) '*2» 
should  be  consulted. 

Within  the  general  topic,  four  separate  tasks  are  specified.  These 
have  been  described  in  detail  in  Reference  1,  but  they  are,  in  brief: 

1.  Investigation  of  fan  and  compressor  design  point  fluid  dynamics, 

2.  Studies  of  compressor  stability  enhancements, 

3.  Fluid  mechanics  of  gas  turbine  operation  in  inlet  flow  distortion 

4.  Investigations  of  three-dimensional  flows  in  highly  loaded 
turbomachines. 

In  addition  to  these  tasks,  the  multi- investigator  contract  also 
encompasses  the  Air  Force  Research  in  Aero  Propulsion  Technology 
program.  The  work  carried  out  in  each  of  the  tasks  will  be  described  in 
the  next  section. 
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TASK  I:  INVESTIGATION  OF  FAN  AND  COMPRESSOR  DESIGN  POINT  FLUID  DYNAMICS 
TASK  IA:  INVERSE  DESIGN  TASK 

During  the  last  three  years,  efforts  on  the  inverse  design  project 
have  been  directed  to  finding  an  acceptable  inverse  method  which  was  appli¬ 
cable  to  the  full  Mach  number  range  encountered  in  high  speed  axial 
compressors.  We  have  found  one  method,  based  on  time  marching  solutions  to 
the  Euler  equations,  which  is  accurate  enough  and  could  potentially  be 
extended  to  three  dimensional  geometries,  see  references  [1.1],  [1.2], 

[1.3],  [1.4].  This  method  operates  in  a  mixed  mode  format  in  which  either 
blade  pressure  distributions  or  geometric  constraints  are  specified.  The 
method  is  unfortunately  far  too  slow,  in  terms  of  computer  time  usage,  for 
practical  design  work  and  in  its  present  form  does  not  properly  adjust 
solid  wall  positions  in  subsonic  flow  regions. 

Work  on  the  time  marching  like  method  is  now  considered  complete  with 
the  publication  of  a  Ph.D.  thesis  on  the  method  by  Tong.  This  thesis 
contains  new  information  about  sources  of  numerical  errors  in  finite 
difference  solutions  and  artificial  smoothing  operators,  a  demonstration 
that  the  simple  inverse  scheme  would  work  in  supersonic  flow  but  not  in 
mixed  subsonic-supersonic  flows,  and  a  survey  and  evaluation  of  existing 
inverse  and  design  methods. 

A  new  direct  or  inverse  method  is  now  being  developed  which  is  many 
times  faster  than  the  old  time  marching  scheme,  is  easily  coupled  to 
boundary  layer  solvers  and  is  more  accurate  than  the  time  marching  style 
solver.  This  new  method  appears  to  retain  all  the  traditional  advantages 
of  streamline  curvature  schemes  while  being  useful  for  subsonic,  supersonic, 


transonic,  and  shocked  flows.  The  scheme  can  be  described  as  a  conserva¬ 
tive  finite  volume  scheme  in  streamline  coordinates  and  has  been  applied  to 
supersonic  duct  flows,  subsonic  and  transonic  duct  flows,  subsonic  cascade 
flows,  and  viscous- inviscid  interactions. 

An  AIAA  paper,  number  84-1643,  was  presented  at  the  June  1984  Fluids 
meeting  in  Snowmass,  CO.  This  paper  presents  theory,  relaxation  methods, 
and  test  examples  for  the  inviscid  applications.  This  paper  is  reproduced 
as  Appendix  1 .  A  paper  on  the  viscous-invicid  interaction  aspects  of  the 
work  has  been  accepted  for  the  Third  Symposium  on  Numerical  and  Physical 
Aspects  of  Aerodynamic  Flows  in  Jan.  1985  at  Long  Beach,  CA.  An  internal 
memo  which  shows  the  present  state  of  the  coupling  work  is  included  as 


Appendix  2 
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CONSERVATIVE  STREAMTUBE  SOLUTION  OF 
STEADY-STATE  EULER  EQUATIONS 

Nark  Drela* 

Massachusetts  Institute  of  Technology 
Cambridge,  Massachusetts 

Michael  Giles* 

Massachusetts  Institute  of  Technology  ■ 
Cambridge,  Massachusetts 

W.  T.  Thompkins,  Jr.® 
Massachusetts  Institute  of  Technology 
Cambridge,  Massachusetts 


Abstract 

This  paper  presents  a  new  method  for  solving 
the  steady  state  Euler  equations.  The  method  is 
similar  to  streamline  curvature  methods  but  has  a 
conservative  finite  volume  formulation  which  ensures 
correct  shock  capturing.  Either  wall  position  or 
wall  pressure  may  be  prescribed  as  boundary  con¬ 
ditions,  permitting  both  direct  and  inverse  calcu¬ 
lations.  In  supersonic  applications  the  solution 
is  obtained  by  space-marching  while  in  subsonic  and 
transonic  applications  iterative  relaxation  methods 
are  used.  Numerical  results  are  given  for: 

a)  Supersonic  diffuser  with  oblique  shocks 
(Direct  calculation) 

b)  Supersonic  jet  entering  still  reservoir 
(Inverse  calculation) 

c)  Subsonic  bump  in  a  channel  with  25%  blockage 
(Direct  and  inverse) 

d)  Subsonic  high-work  turbine  cascade 
(Direct) 

e)  Transonic  bump  in  a  channel  with  12%  blockage 
(Direct  calculation) 

Introduction 

At  present  most  methods  for  the  numerical  cal¬ 
culation  of  steady-state  transonic  solutions  of  the 
Euler  equations  are  based  upon  a  conservative  finite 
volume  approximation  to  the  unsteady  equations  of 
motion  [1].  This  approach  is  conceptually  straight¬ 
forward  since  the  unsteady  equations  are  hyperbolic 
in  time  and  so  the  same  numerical  procedure  can  be 
used  in  regions  where  the  flow  is  locally  subsonic 
or  supersonic.  Using  a  conservative  finite  volume 
formulation  also  guarantees  the  correct  Rankine- 
Hugoniot  shock  jump  relations  regardless  of  the 
details  of  the  shock  calculation. 

The  principal  disadvantage  of  this  approach  is 
that  the  convergence  rate  is  limited  by  the  rela¬ 
tively  slow  propagation  of  w&velike  disturbances 
throughout  the  domain  and  their  reflection  at  the 
boundaries  of  the  computational  domain.  Current 
methods  try  to  overcome  this  limitation  by  a  variety 
of  acceleration  methods  such  as  variable  time  steps, 
implicit  operators  with  larger  time  steps,  and 
multigrid  with  larger  time  steps  on  coarser  grids, 
but  computation  times  remain  an  order-of-magnitude 
larger  than  for  the  solution  of  the  steady-state 
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transonic  full  potential  equation. 

An  older  approach  to  solving  the  subsonic  Euler 
equations  is  the  streamline  curvature  method  which 
can  be  considered  to  be  an  iterative  relaxation 
solution  to  the  equations  of  motion  in  intrinsic 
coordinates  [2,3,4],  In  one  version  [2]  an  initial 
set  of  streamlines  is  guessed  and  the  steady-state 
normal  momentum  equation  is  integrated  in  the  normal 
direction  assuming  known  enthalpy  and  entropy 
variation  along  each  streamtube.  The  resulting 
error  in  the  continuity  equation  then  drives  a  re¬ 
laxation  procedure  to  move  streamlines  towards  their 
correct  positions.  Because  it  is  based  upon  the 
steady  state  equations,  the  relaxation  procedure  in 
the  hyperbolic  supersonic  region  must  be  different 
from  that  in  the  elliptic  subsonic  region. 

One  main  advantage  of  the  streamline  curvature 
method  is  that  it  is  much  faster  than  time-marching 
methods,  and  this  makes  it  much  more  suitable  for 
design  applications  in  which  many  different  geo¬ 
metries  need  to  be  analyzed.  A  secondary  advantage 
is  the  simple  inviscid  treatment  of  adjacent  stream- 
tubes  with  differing  stagnation  properties.  In 
time-marching  calculations  numerical  dissipation 
smears  the  discontinuities. 

The  major  disadvantages  of  the  standard  stream¬ 
line  curvature  method  are  that  the  iterative 
solvers  are  unstable  in  the  supersonic  regions,  and 
they  do  not  conserve  momentum  globally,  so  that  for 
example  the  lift  on  an  airfoil  may  not  equal  the 
change  in  momentum  of  the  fluid.  In  subsonic  cases 
this  momentum  loss  may  be  small,  but  in  transonic 
cases  there  may  be  an  appreciable  momentum  loss  at 
the  shock.  Studies  by  Salas  et  al  [5]  of  solutions 
to  the  full  potential  equation  show  that  momentum 
loss  even  at  relatively  weak  shocks  can  produce 
large  errors  in  shock  position. 

The  method  we  have  developed  is  similar  to  the 
streamline  curvature  method,  but  has  a  conservative 
finite  volume  formulation  which  ensures  the  correct 
treatment  of  shocks.  Like  the  streamline  curvature 
method  it  solves  the  steady-state  equations  and  so 
requires  different  solvers  or  solution  algorithms 
for  purely  supersonic,  purely  subsonic  or  transonic 
problems.  The  supersonic  solver  marches  in  stream- 
wise  steps  downstream  from  a  known  inlet  flow.  The 
subsonic  solver  starts  from  a  guessed  set  of  stream¬ 
line  positions,  solves  each  streamtube  problem  as  a 
separate  quasi  1-0  problem,  and  then  relaxes  all 
the  streamline  positions  to  satisfy  pressure  con¬ 
tinuity  across  streamlines.  The  transonic  solver 
is  similar  to  the  subsonic  solver  but  requires 
numerical  dissipation  in  the  supersonic  region,  and 


(2c) 


relaxes  one  screamline  at  a  time  instead  of  all 
simultaneously.  At  the  present  time  the  supersonic 
and  subsonic  solvers  are  both  robust  and  very  fast, 
but  the  transonic  solver  requires  further  improve¬ 
ment  in  calculation  times. 


Finite  Volume  Formulation  for  Euler  Equations 


The  governing  equations  are  derived  from  the 
integral  form  of  the  mass,  momentum,  and  energy 
conservation  laws. 


q  n  d£  ■  0 
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where  the  integration  is  round  a  closed  curve. 
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These  conservation  equations  are  applied  to 
streamtube  conservation  cells  which  have  the  pro¬ 
perty  that  there  is  no  convection  across  cell  faces 
aligned  with  the  flow  (Figure  1) .  Density  p  and 
speed  q  are  defined  to  be  the  average  values  on 
the  faces  at  the  two  ends  of  the  conservation  cell. 
Also,  the  pressure  in  the  cell  is  assumed  to  be 
linear  and  defined  by  the  average  values  on  the 
four  faces.  For  clarity,  the  average  pressures  on 
the  end  faces  are  denoted  by  p,  and  on  the  stream¬ 
line  faces  by  II. 


Figure  1.  Typical  conservation  cell 


The  flow  direction  at  the  end  faces  is  given  by 
the  unit  vector  s  which  is  defined  by  the  local 
geometry.  The  definition  of  s  in  terms  of  the 
local  cell  nodg  positions  will  be  given  later.  The 
vectors  A  and  B  in  Figure  1  are  the  cell  face  area 
normal  vectors. 

with  the  simplification  that  there  is  no  con¬ 
vection  across  the  streamline  faces,  the  discrete 
conservation  equations  for  mass  and  momentum  are: 

PlVl-Ai  -  02q232*A2  (2a) 

P1A1  +  <0lql®l’Al,ql3l  +  II"B"  “ 

P2A2  +  <02‘l2S2,A2,q2S2  +  (2b) 
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From  the  earlier  assumption  that  the  pressure 
is  linear  over  the  cell,  it  follows  that  the  average 
pressures  on  the  cell  faces  must  be  related  by 

Pl  +  P2  -  +  If  (2d) 

The  assumption  of  no  convection  across  the 
streamline  cell  faces  means  that  the  grid  is  not 
known  a  priori.  The  grid  geometry  must  then  be 
treated  as  an  unknown  which  is  determined  as  part 
of  the  solution.  This  is  a  novel  feature  for  a 
finite  volume  scheme,  although  it  has  long  been 
used  in  streamline  curvature  schemes.  This  feature 
also  makes  the  inverse  and  direc*-  problems  opera¬ 
tionally  equivalent.  In  the  direct  problem,  a 
fixed  wall  position  is  imposed  as  a  boundary  con¬ 
dition  and  the  wall  pressure  is  obtained  as  a  re¬ 
sult.  In  the  inverse  problem,  the  wall  pressure 
is  prescribed  as  the  boundary  condition  and  the 
geometry  is  calculated  as  the  result.  Examples  of 
both  types  of  calculations  will  be  given  later. 


Figure  2.  Geometry  definitions  for  streamline 
curvature  equations 

The  momentum  equation  (2b) ,  which  is  a  vector 
equation,  will  now  be  shown  to  reduce  to  the  set  of 
streamwise  and  normal  momentum  equations  used  by 
the  streamline  curvature  methods,  provided  our 
coordinates,  like  the  s-n  coordinate^,  are^ortho- 
gonal.  More  precisely,  the  vectors  S  and  N  joining 
the  cell  midpoints  in  Figure  2  must  become  perpen¬ 
dicular  and  the  angles  5  and  6  ,  between  ?  and  s, 
and  §2  respectively,  must  vanish  as  the  grid  is 
refined.  taking  the  dot  product  of  (2b)  in  turn 
with  ?  and  N  and  using  (2d) ,  we  obtain  the  following 
streamwise  and  normal  momentum  equations, 

P2q2(  W  (S-’*S)  -  Plql^°l*Al*  <®i*S)  " 

-  jNxS|(p2-Pl>  (3a) 

P2q2 (*2 *A25  (®2,m'  “  piqi**i,Ai> <®i*N>  “ 

-  |nxs|  (if  -  If)  (3b) 

Now  define  m  -  P1'?1<s1*A1)  -  D2q2 (®2‘A2)  • 

With  the  condition  that  ?  and  N  are  orthogonal, 
and  with  the  small  angle  approximations’  sin  6=6 
and  cos  6=1,  equations  (3a)  and  (3b)  reduce  to 


q  -  q  p  —  P 

Since  the  mass  flux  is  constant  along  the  m  2  1  _  _  2  1 

streamtube,  the  enthalpy  conservation  equation  re-  |N|  [sj  Ts| 

duces  to  the  simple  statement  that  total  enthalpy 
is  constant  along  a  streamtube. 


(4a) 
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In  the  Limit  as  the  grid  is  defined,  these 
finite  difference  equations  become  the  streamline 
curvature  equations 


Since  the  streamwise  position  of  the  normal 
grid  lines  is  arbitrary,  some  constraint  on  this 
streamwise  position  must  be  imposed.  At  present, 
the  supersonic  solver  fixes  the  x  positions  for 
simplicity,  leaving  only  the  y's  as  geometric  un¬ 
knowns.  A  more  general  grid  structure  will  be 
treated  in  the  subsonic  solver  section. 


„  aq  dp 

M  r*  *  -  -H- 
03  0S 


0C  0£ 


Thus  our  equations  can  be  interpreted  as  dis¬ 
crete  representations  of  the  streamline  curvature 
equations.  The  essential  point  however  is  that  our 
conservation  form  ensures  that  shocks  will  be 
correctly  captured.  Since  equation  (2d)  was 
necessary  to  obtain  the  above  result,  it  is  inter¬ 
preted  as  a  consistency  condition  on  our  finite 
volume  equations. 

Solution  Procedure  for  Supersonic  Flow 

For  flow  which  is  known  to  be  supersonic 
everywhere,  a  simple  space-marching  solution  pro¬ 
cedure  may  be  used  which  takes  full  advantage  of 
the  hyperbolic  character  of  the  steady  state 
equations.  In  this  solver  the  unknown  geometry 
variables  x,y  are  located  at  the  nodes  of  the  con¬ 
servation  cell  as  shown  in  Figure  3a.  The  s  unit 
vector  giving  the  local  flow  direction  is  defined 
to  be  tangent  to  the  circular  arc  passing  through 
the  three  upstream  midpoints  as  illustrated  in 
Figure  3b.  This  definition  was  chosen  so  that  s 
depends  solely  on  upstream  information,  which  is 
the  correct  domain  of  dependence  for  supersonic 
flow  and  allows  a  space-marching  solution 
procedure. 


At  the  inlet  boundary  all  the  variables  are 
specified.  At  the  top  and  bottom  boundaries  one 
boundary  condition  is  specified.  This  boundary 
condition  can  be  either  a  specified  wall  position 
(direct  case)  or  a  specified  pressure  (inverse 
case) .  with  these  prescribed  boundary  conditions 
the  solver  marches  downstream  from  the  inlet, 
solving  all  of  the  equations  at  each  streamwise 
station  by  using  an  aprpoximate  Newton- Raphson 
linearization  of  the  equations.  This  procedure  is 
similar  to  Keller's  Box  Scheme  [6)  for  the  solution 
of  the  boundary  layer  equations. 

To  reduce  the  number  of  unknowns,  C and  p,,  are 
eliminated  using  the  enthalpy  equation  (2c)  and 
the  pressure  relation  (2d) .  With  these  substitu¬ 
tions  the  equations  for  mass  and  momentum  conser¬ 
vation  are, 


Y q2  s2  (y2  -  y,)  (IT  +TI  -  px> 

R  = - - - 2 - m  -  0 

m  (Y-1XH  -q,/2) 

O  2 

x- momentum 

f*X  =  m(<Jlsl_q2S2  *  +  pl(yl  ”  yi  *  y2  ~  y2)  * 
x  x 

n+(y2“yi)  +  n  iyl'Yp  “  0 


/-momentum 


conservation 

cell 


R  2  mtqjSj^  ~  <I2S2  >  +  -  II+>  (Xj  -  x^)  -  o  (6c) 


In  these  equations  m  is  the  mass  flux  which  is 
constant  for  a  given  streamtube.  s  and  s  are 
the  components  of  the  s  vector.  The  equations  are 
linearized  to  calculate  the  corrections  needed  to 
drive  the  residuals  to  zero. 


Figure  3a. 


Geometry  variable  locations  for 
supersonic  solver 
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Figure  3b.  S  definition  for  supersonic  solver 
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Equations  (7a-c)  for  each  strearatube  cell  are 
assembled  in  a  3x3  block  cridiagonal  form  together 
with  the  two  boundary  conditions!  and  are  solved 
by  standard  methods.  Because  the  equations  are 
fully  linearized,  this  Newton- Raphson  procedure  re¬ 
quires  only  three  or  four  iterations  at  each  stream- 
wise  station  to  converge  machine  accuracy,  p  and 
p.,  can  then  be  calculated  from  the  enthalpy  and 
energy  equations. 


Figure  4  shows  the  final  streamline  grid  and 
Mach  number  contour  plot  for  a  supersonic  channel. 
The  inlet  flow  is  uniform  and  fixed  geometry  is 
specified  on  the  top  and  bottom  walls.  The  agree¬ 
ment  with  the  exact  shock  and  expansion  fan 
positions  is  excellent.  The  shock  has  been 
correctly  captured  without  the  addition  of  any 
artificial  dissipation. 


Figure  4. 


Final  grid  and  Mach  number  contours 
for  supersonic  channel 


Figure  5  shows  the  calculation  of  a  2-D  super¬ 
sonic  jet  entering  a  stagnant  reservoir  filled  with 
a  gas  of  higher  static  pressure.  Before  entering 
the  reservoir  the  flow  is  parallel  with  specified 
top  and  bottom  wall  positions.  After  entering  the 
reservoir  the  top  and  bottom  streamlines  have  con¬ 
stant  specified  pressure.  This  case  physically 
corresponds  to  the  Mach  diamond  pattern  produced 
by  an  overexpanded  supersonic  nozzle  exhausting 
into  a  still  atmosphere.  Both  the  direct  and  the 
inverse  calculations  each  used  only  10  CPU  seconds 
on  a  DEC  PDP  11/70  computer. 


Solution  Procedure  for  Subsonic  Flow 


In  subsonic  flow,  the  steady  state  conservation 
equations  are  elliptic,  since  pressure  signals  can 
propagate  both  upstream  and  downstream.  The  space¬ 
marching  solution  procedure  used  in  the  supersonic 


Figure  5.  Final  grid  and  Mach  number  contours 
for  supersonic  jet 

solver  cannot  properly  represent  the  downstream  in¬ 
fluence  ,  and  hence  an  iterative  relaxation  procedure 
very  similar  to  some  streamline  curvature  solvers, 
is  used  in  the  subsonic  solver.  The  conservation 
cell  and  s  definitions  are  also  different  than  those 
used  in  the  supersonic  solver.  As  shown  in  Figure 
6,  each  conservation  cell  node  is  defined  to  be  at 
the  midpoint  of  a  streamline  segment.  The  s  unit 
vector  lies  along  the  line  connecting  the  midpoints 
of  adjacent  normal  line  segments.  The  boundary  con¬ 
ditions  are  specified  stagnation  pressure,  total 
enthalpy  and  flow  angle  at  the  inlet,  and  specified 
pressure  and  node  position  at  the  exit.  This  last 
condition  is  required  to  determine  the  fraction  of 
the  total  mass  flux  which  is  in  each  streamtube. 


conservation 


Figure  6.  Geometry  variable  locations  and  s 
definition  for  subsonic  solver 

In  the  iterative  relaxation  procedure,  an 
initial  set  of  streamline  positions  is  first  chosen. 
Then  each  streamtube  is  treated  separately  as  a 
quasi-ID  problem  with  equations  (2a),  (2c),  and  (3a) 
being  solved  simultaneously  at  each  streamwise 
station  using  a  Newton-Raphson  technique  incorpora¬ 
ting  the  boundary  conditions.  Equations  (2d)  and 
(3b)  are  then  used  to  calculate  the  pressures  n+, 

II-  on  the  upper  and  lower  surfaces  of  each  stream- 
tube  .  In  the  converged  steady  state  these  must 
match  the  corresponding  pressures  from  neighboring 
streamtubes  so  the  pressure  jump  AH  across  each 
streamline  is  used  to  drive  a  relaxation  of  the 
streamline  positions. 

The  approximate  relationship  between  changes  in 
the  streamline  positions  and  changes  in  the  (1  values 


is  obtained  from  the  linearized  isentropic  pressure- 
area  relation  and  the  pressure-curvature  relation. 


In  these  equations  p,A,n,<  are  the  small  changes 
in  pressure,  cross-sectional  area,  streamline 
pressure  and  curvature,  and  M  is  the  local  Mach 
number.  After  substituting  the  dependence  of  A  and 
<  on  the  changes  in  coordinate  positions  y  we  find 
the  equation  for  the  y  required  to  eliminate  the 
pressure  mismatch  Alt  across  the  streamlines  to  be, 

f(l-M2)u2-52  +  u262]y  - - AH  (10) 

V  n  3  3  nJ  YpM2 

2  2 

Ug  and  5g  are  difference  operators  defined  by 
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2  2 

U  and  Z  are  similarly  defined  operators  in  the 
normal  direction.  Equation  (10)  is  similar  to  the 
streamline  curvature  relaxation  method  used  by 
Keith  et  al  [2] ,  apd  the  iterative  equation  in  the 
solution  of  the  potential  equation  [S].  The 
boundary  conditions  required  for  equation  (10)  are 
y»0  on  the  walls  and  at  the  exit,  and  <5ny-0  at  the 
inlet.  The  latter  corresponds  to  a  fixed  inlet 
flow  angle.  Equation  (10)  is  solved  using  a 
standard  SU3R  technique,  the  streamline  positions 
are  updated  with  y,  and  the  cycle  repeated  until  y 
and  All  become  acceptably  small .  This  SLOR  relax¬ 
ation  technique  limits  the  overall  spectral  radius 
of  the  scheme  to  1-0 (Ax2, Ay2 ) . 


Figure  7a.  Final  grid  and  Mach  number  contours 
for  subsonic  bump 
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Figure  7b.  Relative  stagnation  pressure  loss  for 
subsonic  busp 


A  test  example  is  a  duct  with  a  25%  blockage 
btssp.  The  exit  Mach  number  is  approximately  0.45 
with  a  maximum  Mach  number  in  the  throat  of 
approximately  0.87.  Figure  7a  shows  the  final 
streamline  grid  and  Mach  number  contour  plot.  Be¬ 
cause  the  equations  are  second-order  accurate  and 
no  ntxnerical  dissipation  has  been  needed  to  main¬ 
tain  stability,  the  stagnation  pressure  loss  error 
is  very  small  (Figure  7b) .  The  observed  terminal 
spectral  radius  for  this  40x30  case  was  0.96. 

This  test  case  has  also  been  used  to  check  the 
capabilities  of  the  solver  in  its  inverse  mode. 

The  pressure  on  the  lower  wall  was  set  equal  to 
the  pressure  computed  by  the  direct  method,  and 
the  initial  grid  was  chosen  to  be  the  rectangular 
duct  without  the  bump.  In  approximately  the  same 
number  of  iterations  as  the  direct  case,  the  stream¬ 
lines  and  flow  variables  converged  (to  within 
machine  roundoff)  to  the  values  computed  by  the 
direct  method. 

In  the  case  above,  the  normal  grid  lines  were 
held  at  fixed  x  values  as  in  the  supersonic  case. 
Obviously,  this  is  not  adequate  for  duct  flows 
with  large  turning  angles  or  with  airfoils  with 
blunt  leading  edges.  For  the  turbine  cascade  case 
shown  in  Figure  8  a  more  suitable  constraint  on 
the  normal  grid  lines  was  used.  During  the  stream¬ 
line  relaxation  step,  the  nodes  are  moved  parallel 


to  the  local  normal  lines.  Then  the  nodes  are 
moved  along  the  streamlines  by  an  algorithm  which 
minimizes  the  local  shearing  of  the  grid  and  vari¬ 
ations  in  streamwise  resolution.  The  Kutta  condition 
is  imposed  by  shifting  the  entire  exit  plane  ver¬ 
tically  until  the  values  of  II  on  either  side  of  the 
trailing  edge  are  equal. 

All  the  cases  shown  above  have  either  hardwall 
(prescribed  geometry)  or  freewall  (prescribed 
pressure)  boundary  conditions,  with  our  scheme,  it 
is  also  possible  to  prescribe  a  "softwall"'  boundary 
condition  where  the  wall  is  allowed  to  deform  under 
the  effect  of  the  wall  pressure  according  to  some 
prescribed  relationship.  Because  the  displacement 
surface  of  a  boundary  layer  can  be  described  in  this 
manner,  the  softwall  boundary  condition  makes 
efficient  and  straightforward  boundary  layer  coupling 
possible  even  in  strongly  coupled  cases  with  separa¬ 
tion  bubbles.  The  details  and  numerical  examples  of 
the  coupling  between  this  Euler  solver  and  a  com¬ 
pressible  boundary  layer  code  (7)  will  be  presented 
in  a  future  paper. 

Transonic  Solver 

The  transonic  solver  is  still  under  development. 
Currently  it  is  not  fast  as  the  subsonic  solver,  but 
futvre  planned  changes,  including  multigrid,  should 
greatly  increase  its  speed.  The  conservation  cell 
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U  >  max 


The  iterative  solution  procedure  is  to  guess  an 
initial  streamline  grid  and  then  sweep  up  and  down 
in  the  normal  direction,  relaxing  each  streamline 
by  assuming  neighborino  streamlines  are  fixed  and 
so  solving  the  steady  state  equations  for  the  two 
streamtubes  on  either  side  of  the  streamline  using 
a  Newton- Raphson  procedure.  Figure  9a  shows  the 
final  grid  and  Mach  number  contour  plot  for  a 
partially  choked  duct  with  a  12%  blockage  bump. 
Figurt  j.  gives  the  relative  stagnation  pressure 
losses. 


Figure  9a.  Final  grid  and  Mach  number  contours 
for  transonic  bump. 


Figure  8. 


Final  grid  and  Mach  number  contours 
for  turbine  case 


and  s  definitions  are  the  same  as  for  the  subsonic 
solver.  To  capture  shocks,  however,  "artificial 
compressibility"  is  introduced  in  the  mass 
continuity  equation. 

^1^1  ^*1**1  *  *  °2<J2^2*A2) 
where 

s2  “  (i-u)e2 

and  o  is  similarly  defined.  This  artificial  com¬ 
pressibility  was  first  used  by  A.  Eberle  [8]  in 
solving  the  transonic  full  potential  equation,  and 
has  been  used  more  recently  by  vtomom  (91  in 
solving  the  quasi  1-D  Euler  equations,  u  is  a 
function  of  the  Mach  number  M  and  to  ensure  numeri¬ 
cal  stability  must  satisfy 


Figure  9b.  Relative  stagnation  pressure  loss 
for  transonic  bump 

Conclusions 

This  paper  has  presented  a  new  method  for 
solving  the  steady-state  Euler  equations.  The 
numerical  examples  presented  show  its  capability  to 
handle  supersonic,  subsonic  and  transonic  flows  with 
direct  or  inverse  boundary  conditions.  In  addition 
the  cascade  example  shows  the  ability  to  determine 
the  correct  circulation  around  an  airfoil  by  ex¬ 
plicitly  applying  a  Kutta  condition  at  the  trailing 
edge.  Further  work  is  required  to  improve  the 
accuracy  and  speed  of  the  transonic  solver,  and  to 
formulate  a  transonic  inverse  procedure. 
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Abstract 


This  report  presents  preliminary  results  obtained  using  a  new 
method  for  solving  the  steady  state  Euler  equations.  The  method  is 
similar  to  streamline  curvature  methods  but  has  a  conservative  finite 
difference  formulation  which  ensures  correct  shock  capturing.  In 
supersonic  applications  the  solution  is  obtained  by  space-marching  while 
in  subsonic  applications  an  iterative  elliptic  relaxation  procedure, 
similar  to  potential  methods,  is  used.  Numerical  results  are  given  for; 

a)  Supersonic  jet  entering  still  reservoir  (Inverse  calculation) 

b)  Subsonic  Ni's  bump  with  40%  blockage  (Direct  and  inverse) 

c)  Subsonic  bump  with  boundary  layer  coupling  to  calculate  a  separation 
bubble 

d)  Transonic  Ni's  bump  with  10%  blockage  (Direct  calculation) 


14 


1 )  Introduction 


At  present  most  methods  for  the  numerical  calculation  of  steady- 
state  transonic  solutions  of  the  Euler  equations  are  based  upon  a 
conservative  finite  volume  approximation  to  the  unsteady  equations  of 
motion?.  This  approach  has  two  good  features.  The  first  is  that  it 
is  conceptually  straightforward  since  the  unsteady  equations  are 
hyperbolic  in  time  and  so  the  same  numerical  procedure  can  be  used  in 
both  the  subsonic  and  supersonic  regions.  The  second  is  that  the 
conservative  finite  volume  formulation  guarantees  the  correct  Rankine- 
Hugoniot  shock  jump  relations  regardless  of  the  details  of  the  shock 
calculation.  The  principal  disadvantage  is  that  the  convergence  rate  is 
limited  by  the  relatively  slow  propagation  of  vavelike  disturbances 
throughout  the  domain  and  their  reflection  at  the  boundaries  of  the 
computational  domain.  Current  methods  try  to  overcome  this  limitation 
by  a  variety  of  acceleration  methods  such  as  variable  time  stepB, 
implicit  operators  with  larger  time  steps,  and  multigrid  with  larger 
time  steps  on  coarser  grids,  but  computation  times  remain  an 
order-of-magnitude  larger  than  for  the  solution  of  the  steady-state 
transonic  full  potential  equation,  using  iterative  relaxation  methods, 
enhanced  by  multigrid. 

An  older  approach  is  the  Streamline  Curvature  method  which  is  an 
iterative  relaxation  procedure  which  solves  the  equations  of  motion  in 
intrinsic  coordinates1 >5,6.  in  one  version1  an  initial  set  of 
streamlines  is  guessed,  the  steady-state  normal  momentum  equation  is 


integrated  in  the  normal  direction  aasuming  known  enthalpy  and  entropy 
variation  along  each  streamtube.  The  resulting  error  in  the  continuity 
equation  then  drives  a  relaxation  procedure  to  move  streamlines  towards 
their  correct  positions.  Because  it  is  based  upon  the  steady  state 
equations,  the  relaxation  procedure  in  the  hyperbolic  supersonic  region 
is  different  from  that,  in  the  elliptic  subsonic  region.  The  main  advan¬ 
tage  of  the  streamline  curvature  method  is  that  it  is  much  faster  than 
time-marching  methods,  and  this  makes  it  much  more  suitable  for 
industrial  design  applications  in  which  many  different  geometries  need 
to  be  analyzed.  A  secondary  advantage  is  the  simple  inviscid  treatment 
of  adjacent  streamtubes  with  differing  stagnation  properties.  In  time- 
marching  calculations  numerical  dissipation  smears  the  discontinuities. 
The  disadvantages  of  the  standard  streamline  relaxation  method  are  that 
the  iterative  solvers  have  problems  with  numerical  stability  in  the 
supersonic  regions,  and  they  do  not  conserve  momentum  globally,  so  that 
for  example  the  lift  on  an  airfoil  may  not  equal  the  change  in  momentum 
of  the  fluid.  In  subsonic  cases  this  momentum  loss  may  be  small,  but  in 
transonic  cases  there  may  be  an  appreciable  momentum  loss  at  the  shock. 
Studies  by  Salas  et  al4  of  solutions  to  the  full  potential  equation 
show  that  momentum  loss  even  at  relatively  weak  shocks  can  produce  large 
errors  in  shock  position. 

The  method  we  have  developed  is  similar  to  the  streamline  curvature 
method,  but  has  a  conservative  finite  volume  formulation  which  ensures 
the  correct  treatment  of  shocks.  Like  the  streamline  curvature  method 
it  solves  the  steady-state  equations  and  so  requires  different  solvers 


or  solution  algorithms,  for  pursly  supsrsonic,  pursly  subsonic,  or 


transonic  problems.  The  supersonic  solver  marches  in  streamwise  steps 
downstream  from  a  known  inlet  flow.  The  subsonic  solver  starts  from  a 
guessed  set  of  streamline  positions,  solves  each  streamtube  problem  as  a 
separate  quasi  1 -D  problem,  and  then  relaxes  all  the  streamline 
positions  to  satisfy  pressure  continuity  across  streamlines.  The  tran¬ 
sonic  solver  is  similar  to  the  subsonic  solver  but  requires  numerical 
dissipation  in  the  supersonic  region,  and  relaxes  one  streamline  at  a 
time  instead  of  all  simultaneously.  At  the  present  time  the  supersonic 
and  subsonic  solvers  are  both  robust  and  very  fast,  but  the  transonic 
solver  requires  further  improvement. 

2)  Supersonic  Solver 

The  geometry  and  location  of  flow  variables  for  a  typical  "cell"  of 
a  streamtube  is  shown  in  Figure  1 .  The  geometry  is  defined  at  the 
comers  of  each  cell.  Density  p,  speed  q,  and  pressure  p  are  defined  at 
the  midpoints  of  the  two  ends  of  the  cell  and  the  auxiliary  pressure  n 
is  defined  on  the  streamline  faces  of  the  cell,  s  is  a  unit  vector  in 
the  direction  of  the  velocity,  and  is  defined  to  be  tangent  to  the 
circular  arc  passing  through  the  three  upstream  midpoints  as  illustrated 
in  Figure  2.  At ,  A,,  B~ ,  B"  are  area  vectors  normal  to  the  cell  faces. 

The  finite  volume  formulation  of  the  steady  equations  is 
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In  addition  there  is  an  auxiliary  pressure  equation  which 
constrains  the  average  of  the  pressures  on  the  two  pairs  of  opposing 
faces  to  be  equal. 

Pi  +  Pi  *  n+  +  n  (id) 

At  the  inlet  boundary  all  the  variables  need  to  be  specified.  At 
the  top  and  bottom  boundaries  one  boundary  condition  needs  to  be 
specified,  and  at  present  the  code  accepts  either  a  specified  wall 
position  (direct  case)  or  a  specified  pressure  (inverse  case).  With 
these  prescribed  boundary  conditions  the  solver  marches  downstream  from 
the  inlet,  solving  all  of  the  equations  at  each  streamwise  station  by 
using  an  approximate  Newton-P.aphsen  linearization  of  the  equations, 
which  requires  typically  about  four  iterations  at  each  station.  The 
process  is  similar  to  Keller's  box  method  for  the  solution  of  the 
boundary  layer  equations. 

Figure  3  shows  a  Mach  number  contour  plot,  and  the  final  streamline 
grid  for  a  2-0  supersonic  jet  entering  a  stagnant  reservoir  filled  with 
a  gas  of  higher  static  pressure.  Before  entering  the  reservoir  the  flow 
is  parallel  with  specified  top  and  bottom  walls.  After  entering  the 
reservoir  the  top  and  bottom  streamlines  have  constant  specified 


pressure .  This  case  physically  corresponds  to  the  Mach  diamond  pattern 
produced  by  a  supersonic  nozzle  exhausting  into  a  still  atmosphere  and 
took  only  10s  on  a  PDP11/70  to  achieve  convergence  to  machine  roundoff. 

3)  Subsonic  Solver 

Initially  our  subsonic  solver  used  the  same  geometry  and  location 
of  flow  variables  as  the  supersonic  solver.  Unfortunately  there  were 
problems  with  numerical  stability  and  unacceptably  high  errors.  These 
problems  were  overcome  by  going  to  the  slightly  different  geometry  shown 
in  Figure  4,  without  needing  to  resort  to  numerical  dissipation.  There 
are  now  two  grid  structures,  a  relaxation  grid  with  coordinates  defined 

at  its  nodes,  and  a  conservation  grid  with  its  comer  nodes  defined  to 
be  the  midpoint  of  the  corresponding  streamline  faces  on  the  relaxation 

grid.  The  s  unit  vector  in  the  direction  of  the  velocity  is  defined 
to  be  parallel  to  a  line  connecting  midpoint  of  normal  faces  on  the 
relaxation  grid  as  shown  in  Figure  5.  The  steady  state  conservation 
equations  are  exactly  the  same  as  for  the  supersonic  solver  but  the 
solution  procedure  is  totally  different.  Since  the  equations  are 
elliptic,  an  iterative  relaxation  procedure,  very  similar  to  some 
streamtube  curvature  solvers,  is  required.  First  an  initial  set  of 
streamline  positions  is  chosen.  Then  each  streamtube  is  treated 
separately  and  Equations  (la-id)  are  solved  simultaneously  at  each 
streamwise  station  using  a  Newton-Raphson  technique  incorporating  the 
boundary  conditions,  specified  stagnation  pressure  total  enthalpy  and 
exit  pressure.  These  calculations  produce  pressures  n  +  ,  H”  on  the 


upper  and  lower  surfaces  of  each  streamtube.  In  the  converged  steady 
state  these  oust  match  the  corresponding  pressures  from  neighboring 
streamtubes  so  the  pressure  jump  AH  across  the  streamline  drives  a 
relaxation  of  the  streamline  positions.  Assuming  locally  isentropic 
behavior,  influence  coefficients  are  calculated  for  the  effect  of 
geometry  changes  on  the  local  streamtube  areas  and  curvature,  and  hence 
on  the  II  values.  The  resultant  relaxation  equation  for  the  case  in 
which  the  "normal"  lines  are  orthogonal  is 

( 1-M*  )6*  +  6*  y  - - r-  AH  (2) 

8  n  pu  A 

In  this  equation  M  is  the  Mach  number,  y  is  the  normal  correction 
in  the  streamline  position,  A  is  the  streamtube  area  and  ah  is  the 
forcing  pressure  mismatch  across  the  streamline.  6s  and  6n  are 
finite  difference  approximations  to  second  derivatives  in  the  streamvise 
and  normal  directions.  Thus  this  equation  is  very  similar  to  both  the 
transonic  small  disturbance  equation  and  the  streamtube  curvature 
relaxation  method  used  by  at  al  • .  We  solve  (2)  using  a  standard 

SLOR  technique,  update  the  streamline  positions  and  then  repeat  the 
cycle  until  y  and  ah  become  acceptably  small. 

The  test  example  is  a  duct  with  a  40%  blockage  bump.  The  inlet 
Mach  number  is  approximately  0.27  with  a  maximum  Mach  number  in  the 
throat  of  approximately  0.52.  Figure  6  shows  the  initial  and  final 
streamline  grids.  To  increase  the  resolution  at  the  stagnation  points 
two  narrow  streamtubes  are  used  next  to  the  lower  wall  and  there  are 


additional  normal  lines  at  the  front  and  back  of  the  bump.  The  only 
discemable  difference  between  the  initial  and  final  grid  is  a 
thickening  of  the  streamtubes  near  the  stagnation  point.  Contour  plots 
of  the  Mach  number  and  relative  stagnation  pressure  loss  are  shown  in 
Figure  7.  Because  the  equations  are  second-order  accurate  and  no 
numerical  dissipation  has  been  needed  to  maintain  stability,,  the 
stagnation  pressure  loss  is  ouch  lower  than  for  time-marching  codes2. 

This  test  case  has  also  been  used  to  check  the  capabilities  of  the 
solver  in  its  inverse  mode.  The  pressure  on  the  lower  wall  was  set 
equal  to  the  pressure  computed  by  the  direct  method,  and  the  initial 
grid  was  chosen  to  be  the  rectangular  duct  without  the  bump.  In  approx¬ 
imately  50  iterations  the  streamlines  and  flow  variables  converged  (to 
within  machine  roundoff)  to  the  values  computed  by  the  direct  method. 

The  ability  to  prescribe  wall  position  and/or  wall  pressure  in  our 
formulation  makes  efficient  boundary  layer  (BL)  coupling  possible  even 
in  strongly  coupled  cases  with  separation  bubbles  (Figures  8,  9a, b). 

For  this  example,  a  finite  difference  BL  code  using  a  Shifted  Box 
Scheme8  was  used.  Simultaneous  BL  calculations  and  downstream  relax¬ 
ation  sweeps  are  used  to  drive  the  Euler  and  BL  edge  velocities  together 
at  the  displacement  surface.  During  the  relaxation  of  the  inviscid 
flow,  a  hybrid  boundary  condition  prescribing  a  linear  combination  of 
edge  pressure  and  position  simultaneously  drives  these  quantities 
towards  the  matching  condition.  The  behavior  of  the  BL  is  locally 
modeled  by  a  u,  -  {•  influence  coefficient  which  is  a  byproduct  of  the 
BL  calculation. 


The  coupled  solution  in  Figures  9a, b  required  120  relaxation  sweeps 
and  12  Bi>  calculations  to  fully  converge.  The  purely  inviscid  case  in 
Figure  8  required  50  relaxation  sweeps.  This  preliminary  coupling  case 
indicates  that  coupled  calculations  will  not  require  drastically  more 
computation  time  than  purely  inviscid  calculations. 

4)  Transonic  Solver 

The  transonic  solver  is  still  under  development.  The  current 
solver  is  not  very  fast  or  robust  but  is  sufficient  to  demonstrate  the 
feasibility  of  this  approach  for  transonic  cases.  The  grid  geometry  is 
the  same  as  the  supersonic  solver.  At  supersonic  points  the  supersonic 
i  definition  is  used.  At  subsonic  points  s  is  defined  to  be  tangent 
to  the  circle  passing  through  the  neighboring  midpoints  as  shown  in 
Figure  1 .  The  main  feature  of  the  transonic  steady  state  equations  is 
the  introduction  of  an  "artificial  compressibility"  term  in  the 
continuity  equation.  The  mass  conservation  equation  is 

(81*N  }  > 

where 

p,  ■  ( 1-u )p,  +  u  P, 
and  pL  is  similarly  defined. 

This  artificial  compressibility  was  first  used  by  A.  Eberle7  in 
solving  the  transonic  full  potential  equation,  and  has  been  used  more 
recently  by  Womom  [3]  in  solving  the  quasi  1-D  Euler  equations,  u  is 
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a  function  of  tha  Mach  number  M  and  to  ensure  numerical  stability  must 
satisfy 


u  >  max 


The  iterative  solution  procedure  is  to  guess  an  initial  streamline 
grid  and  then  sweep  from  side  to  side  in  the  normal  direction,  relaxing 
each  streamline  by  assuming  neighboring  streamlines  are  fixed  and  so 
solving  the  steady  state  equations  for  the  two  streamtubes  on  either 
side  of  the  streamline  using  a  Newton-Raphson  procedure. 

Figure  10  shows  the  final  grid  and  Mach  number  contour  plot  for  a 
fully  choked  duct  with  a  25%  blockage  bump. 
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FIGURE  4 .  Finite  Volume  Geometry  for  Subsonic  Solver 


5.  Definition  of  Subsonic  Unit  Velocity  Vector  S. 
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FIGURE  8.  Mach  number  contour  plot  without  boundary  layer  interaction 


0.  350 
0.  325 
0.  300 
0.275 
0.250 
0.225 


0.400 


contour 


TASK  IB:  LOSS  MECHANISMS  AND  LOSS  MIGRATION  IN  TRANSONIC  COMPRESSORS 
Introduction 

The  objective  of  this  study  is  to  quantify  the  understanding  of  loss 
production  versus  loss  migration  by  detailed  time-resolved  radial 
streamline/pathline  measurements  in  the  AFAPL  High  Through  Flow 
Compressor,  a  transonic  compressor  stage,  in  the  MIT  Blowdown  Compressor. 

Streamline  tracing  in  a  transonic,  unsteady  compressor  flow  requires 
that  streamlines  or  pathlines  be  identified  in  the  presence  of 
simultaneous  first  order  nonperiodic  changes  in  total  pressure  and 
temperature . 

The  basic  technique  is  to  inject  a  thin  ring  of  the  tracer  gas  upstream 
of  the  compressor  rotor  and  then  survey  with  a  concentration  probe  at  the 
rotor  exit  to  map  the  radial  transport  of  the  tracer  gas.  The  injection 
ring  is  then  moved  radially  and  the  experiment  repeated.  In  this  manner 
the  time-dependent  radial  transport  within  a  transonic  compressor  rotor  can 
be  mapped  simultaneously  with  a  timeresolved  measurement  of  rotor  loss 
(from  a  cotemperal  and  pressure  measurement). 

The  dual  wire  aspirating  probe  will  be  used  for  the  simultaneous 
measurement  of  time-resolved  concentration,  temperature,  and  pressure 
possible  at  frequencies  up  to  20  KHz.  The  probe  consists  of  two  constant 
temperature  hot  wires  operated  at  different  overheat  ratios  in  a  channel 
with  a  choked  exit.  Thus  the  flow  past  the  wires  is  at  constant  Mach 
number.  The  heat  flux  (and  thus  the  electrical  signal)  from  the  wires  is  a 
function  only  of  freestream  total  pressure  and  temperature  and  gas  composi¬ 
tion.  Previously,  the  probe  was  used  in  a  flow  of  constant  composition. 


Thus,  the  signals  from  the  two  wires  were  sufficient  to  permit  calculation 


of  the  two  unknowns,  total  pressure  and  temperature.  A  piggyback  total 
pressure  sensor  was  used  as  a  redundant  check.  In  this  investigation,  the 
pressure  sensor  is  used  to  provide  the  third  signal  to  allow  the  three 
independent  variables  (concentration,  temperature,  and  pressure)  to  be 
derived  at  each  instant  in  time. 

Progress  to  Date 

Work  to  date  has  concentrated  on  tracer  gas  selection,  probe 
calibration,  and  injector  design.  The  selection  of  the  proper  tracer  gas 
is  important.  Ideally,  the  tracer  should  have  the  same  density  as  the 
main  flow  while  differing  in  those  fluid  properties  important  to  the 
probe  operation:  thermal  conductivity,  viscosity,  ratio  of  specific 
heats,  and  molecular  weight.  We  have  examined  the  suitability  of  several 
gases  and  gas  mixtures  in  regard  to  the  above  criteria  and  concluded  that 
the  two  most  promising  for  use  in  the  Argon-Freon- 1 2  atmosphere  of  the  MIT 
Blowdown  Compressor  are  carbon  dioxide  and  Freon  153A.  Krypton  is  also 
suitable  but  quite  expensive.  Figure  1  illustrates  the  calculated  voltage 
output  from  the  aspirating  probe  as  a  function  of  gas  temperature  and 
CO2  fractional  concentration  in  the  Argon-Freon  mixture  used  in  the  MIT 
Blowdown  Compressor  Facility. 

It  is  clear  from  the  analysis  that  the  success  of  the  experiment 
hinges  on  the  overall  stability  of  the  aspirating  probe  since  the  variations 
due  to  concentration  changes  are  relatively  small,  representing  the 
difference  between  two  large  numbers.  Preliminary  laboratory  calibration 
tests  with  CO2  look  promising.  These  tests  are  continuing,  to  both 


establish  the  calibration  and  data  reduction  procedures  and  to  assess  the 
long  term  probe  stability. 

The  gas  injector  has  been  designed  and  is  now  being  fabricated.  The 
design  goal  is  the  injection  of  a  1  mm  thick  sheet  of  tracer  gas  with  the 
same  velocity  and  static  pressure  as  the  main  flow  in  order  to  minimize  the 
turbulent  mixing  of  this  sheet  as  gas  propagates  downstream.  The  injector 
(figure  2)  consists  of  a  NACA  654-021  airfoil  limited  to  two  blade 
passages  in  circumferential  extent  in  order  to  reduce  the  amount  of 
injection  gas  required.  The  hollow  airfoil  contains  a  plenum  (instrumented 
with  a  pressure  transducer)  feeding  a  0.8  mm  wide  slot  at  the  trailing 
edge.  The  injector  is  supported  by  two  radial  struts  of  similar 
cross-section  which  carry  the  injection  gas.  The  entire  assembly  is 
mounted  through  the  tip  casing  which  permits  the  radial  position  to  be 
varied.  O-rings  provide  the  seal  at  this  point.  The  vibrational 
frequencies  of  the  injector  assembly  has  been  calculated  to  insure  that  no 
resonances  lie  within  the  range  that  can  be  forced  by  the  rotor. 

An  attempt  was  made  to  calculate  the  extent  of  the  radial  mixing  due 
to  turbulent  shear  as  the  injectant  sheet  is  convected  downstream.  Because 
this  is  such  a  strong  function  of  the  precise  pressure  and  velocity  match 
at  the  Injector  trailing  edge  and  because  various  models  give  such 
disparate  results,  the  only  believable  figures  will  be  those  from  measure¬ 
ments.  These  will  be  done  by  surveying  across  the  compressor  annulus  at 
the  rotor  exit  station  with  the  blade  rows  removed,  thus  establishing  a 
baseline  for  comparison  with  the  rotor  exit  measurements. 

It  is  anticipated  that  most  of  the  calibration  work  will  be  completed 
by  the  time  the  injector  assembly  is  fabricated.  The  baseline  measurements 
described  above  will  then  begin. 
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Aspirating  probe  output  voltage  as  a  function  of  C02  concentration  at 
different  freestream  total  temperatures.  The  total  pressure  is  0.3  atm 


TASK  II:  COMPRESSOR  STABILITY  ENHANCEMENT  USING  HUB/CASING  TREATMENT 


The  overall  goal  of  this  study  is  the  understanding  of  the  basic 
mechanism(s)  by  which  a  grooved  casing  or  hub  is  able  to  significantly 
extend  the  stable  flow  range  of  a  turbomachine.  A  key  concept  in  this  is, 
in  our  view,  a  quantitative  appreciation  for  the  fluid  mechanics  of  the 
interchange  of  momentum  between  the  casing  and  the  flow  in  the  endwall 
region  of  the  blade  passage. 

In  our  first  series  of  experiments,  carried  out  over  the  past  several 
years,  we  showed  that  a  moving  grooved  hub  could  have  a  very  potent  effect 
on  stator  stall  margin  and  pressure  rise  [2.1].  (This  was  the  first  time 
this  had  been  demonstrated. )  This  result  meant  that  we  could  examine  the 
relevant  phenomena  using  the  stator  passage/moving  hub  configuration.  The 
present  series  of  experiments  are  aimed  much  more  at  the  basic  phenomena 
associated  with  hub  and  casing  treatment,  and  involve  a  detailed  examina¬ 
tion  of  the  three-dimensional,  unsteady  flow  within  the  blade  passage. 

This  work  will  form  part  of  the  M.S.  Thesis  of  Mr.  Mark  Johnson,  who  is  in 
the  AFRAPT  Program. 

In  previous  progress  reports,  we  mentioned  the  design,  development, 
and  fabrication  of  the  instrumentation  and  the  necessary  software  to  enable 
us  to  examine  in  detail  the  unsteady  flow  in  the  stator  passage.  He  also 
mentioned  the  design  of  an  optical  detector  system  to  interact  with  our  LSI 
11/23  microcomputer  to  monitor  the  relative  position  of  the  grooves 
accurately  during  the  experiment.  In  the  6-month  period  following  that 
progress  report,  we  have  essentially  completed  the  tasks  mentioned  above. 

In  addition,  we  have  also  accomplished  the  following  tasks: 


( 1 )  With  the  newly  designed  and  built  inlet  contraction  section  mounted 
on,  measurements  of  the  inlet  velocity  profile  upstream  of  the  single  stage 
low  speed  compressor  were  taken  to  compare  with  the  previous  results  as  a 
check.  No  substantial  difference  in  the  results  was  observed  except 

that,  as  might  be  expected  because  of  the  longer  flow  path,  the  boundary 
layer  on  the  casing  is  thicker  than  before  (i.e.,  with  the  old  inlet 
section  on ) . 

(2)  The  compressor  characteristics,  pressure  rise  versus  flow,  were 
retaken  both  for  the  cases  of  smooth  hub  and  treated  hub  to  compare  with 
previous  results  as  a  check.  This  data  is  shown  as  Figure  2.1.  The  large 
differences  between  the  smooth  wall  and  the  grooved  wall  pressure  rise  as 
well  as  stall  onset  points  (essentially  at  the  peaks  of  the  curves)  are 
evident.  Although  the  relative  positions  of  the  curves  are  close  to  those 
measured  previously,  there  is  a  difference  (of  about  5%)  in  the  absolute 
levels  of  flow  and  pressure  from  previous  measurement  for  both  smooth  wall 
and  hub  treatment.  This  may  be  due  to  the  larger  inlet  boundary  layer 
resulting  in  a  changed  profile  presented  to  the  stator.  However,  because 
of  the  departure  of  Mr.  Johnson  for  Garrett  (where  he  is  working  for  the 
summer,  as  part  of  the  AFRAPT  program)  we  have  not  resolved  this  with  any 
degree  of  precision. 

(3)  Preliminary  experiments  were  carried  out  to  validate  and  test  the 
overall  data  acquisition  system.  This  includes  the  traverse  system  for 
probe  (hot  wire  or  pressure)  with  three  degrees  of  freedom  (linear 
translation,  angular  motion  in  a  radial  axial  plane  in  the  stator  blade 
passage,  and  rotation  on  its  own  axis)  and  the  optical  detector  system,  for 


monitoring  the  groove  position. 

An  inclined  hot-vire  probe  was  used  to  obtain  high  response  data  at  a 
particular  point  in  the  blade  passage,  from  which  the  three  velocity 
components  at  that  point  can  be  reduced.  This  probe  had  been  separately 
calibrated  in  a  GTL  jet  facility.  The  spatial  resolution  achievable  with 
this  current  hot  wire  probe  is  2  to  3  mm. 

The  data  matrix  grid  can  be  described  by  a  three-dimensional 
co-ordinate  system  (I, J,K)  fixed  to  the  laboratory  frame.  Bach  1  refers  to 
a  plane  concentric  with  the  hub  surface  (for  instance  1-1  refers  to  a 
plane  2.25  mm  from  the  hub,  I  ■  2  to  that  12.25  mm  from  the  hub,  etc.). 

Each  J  refers  to  a  plane  parallel  to  a  mean-chord  plane  with  J  -  1  plane 
being  closest  to  the  blade  suction  surface.  Finally,  each  K  refers  to  a 
plane  parallel  to  the  leading  end  trailing  edge  plane  with  K  -  1  plane 
being  the  leading  edge  plane. 

For  this  set  of  preliminary  experiments,  the  data  matrix  chosen  is 
such  that  on  each  I  plane  (I  varies  from  1  to  5,  I  -  5  is  roughly  the 
mid-span  plane)  there  is  a  4  by  7  grid  (i.e.,  J  «  1  to  4  while  K  -  1  to  7). 
A  separate  set  of  data  was  taken  at  grid  points  outside  the  above  data 
matrix  grid  for  ascertaining  the  periodicity  of  the  flow  from  passage  to 
passage.  In  the  preliminary  experiments,  hot  wire  data  were  acquired  for 
runs  with  smooth  hub  wall  as  well  as  groove- treated  hub  wall,  with  the 
compressor  operating  at  2600  rpm  and  at  flow  coefficients  varying  from  0.24 
to  0.37.  The  hot  wire  data  have  not  yet  been  reduced  to  yield  the  three- 
dimensional  flow  field  in  the  blade  passage  region  close  to  the  smooth  and 
treated  hub  surfaces.  However,  it  does  appear  from  this  set  of  experi- 
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Bents,  that  the  overall  data  acquisition  system  performs  satisfactorily. 

We  do  not  as  yet  rule  out  any  future  modifications  and  refinements  to  the 
data-acquisition  system  if  ve  find  it  necessary. to  do  so. 

We  hope  that,  in  the  next  few  months,  efforts  will  be  expended  to 
obtain  the  three-dimensional  velocity  field  in  the  passage  from  this  hot¬ 
wire  data,  and  to  examine  possible  three-dimensional  presentations  of  this 
information.  Such  a  presentation  would  enhance  our  ability  to  view  the 
results  in  a  more  interpretative  and  constructive  way.  This  preliminary 
set  of  experiments  should  also  enable  us  to  refine  and  to  plan  our 
measurements  of  the  flow  in  the  next  set  of  experiments.  When  this  set  of 
measurements  has  been  carried,  out,  we  also  hope  to  take  up  again  our 
theoretical  (modelling)  work  on  this  problem.  This  is  an  area  which  we 
have  not  had  time  to  pursue  over  the  past  year  but  one  that  is  a  crucial 
part  of  a  balanced  research  effort  on  this  topic. 

Reference 
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TASK  III:  INLET  VORTEX  FLOW  DISTORTIONS  IN  GAS  TURBINE  ENGINES 

During  the  past  six  months,  we  have  nearly  completed  the  experimental 
part  of  our  investigations  of  this  phenomenon  (at  least  the  first  phase  - 
the  vortex  formation  process).  These  have  provided  not  only  quantitative 
information  on  vortex  strength,  but  also  insight  into  the  basic  vortex 
structure  associated  with  this  type  of  flow.  A  paper  on  this  work  is 
currently  being  prepared  for  submission .  A  draft  of  this  paper  is 


presented  below 


INTRODUCTION 


Gas  turbine  engine  operation  near  the  ground  is  often  associated  with 
the  presence  of  a  strong  vortex  which  stretches  between  the  ground  and  the 
engine  inlet.  This  "inlet  vortex"  (or  ground  vortex)  can  be  of  concern 
because  of  the  effects  it  can  have  on  the  engine.  These  include  foreign 
object  damage,  increased  blade  erosion  due  to  dirt  ingestion  [1],  and 
compressor  surge  [2]. 

Investigations  of  the  inlet  vortex  have  identified  two  basic  mechanisms 
that  are  responsible  for  its  formation  [3].  The  first  is  the  amplification 
of  the  vertical  component  of  ambient  vorticity,  i.e.,  vorticity  created  far 
us tr earn  of  the  inlet,  as  the  vortex  lines  are  convected  into  the  inlet.  A 
detailed  description  of  this  process,  including  calculations  of  the  drift 
of  these  vertical  vortex  lines  in  the  three-dimensional  flow  produced  by  an 
inlet  near  a  ground  plane ,  is  given  in  [ 3 } . 

The  second  mechanism  for  vortex  formation,  however,  does  not  depend  on 
the  presence  of  ambient  vorticity.  This  process  can  occur  in  a  flow  that 
is  irrotational  far  upstream,  for  an  inlet  in  crosswind.  In  this  situation, 
because  of  the  variation  in  circulation  along  the  inlet,  a  trailing  vortex 
(from  the  rear  aide  of  the  inlet)  is  also  present.  It  is  this  second 
mechanism,  which  appears  to  be  less  well  understood,  that  forms  the  focus 
of  the  present  paper. 

In  previous  investigations,  the  flow  regimes  encountered  with  an 
inlet/ground  plane  configuration  were  found  to  depend  strongly  on  two 
parameters:  1)  the  ratio  of  inlet  (average)  velocity  to  upstream  velocity, 
Ui/U.,  and  2)  the  ratio  of  inlet  centerline  height  to  inlet  inside 
diameter,  H/D.  For  an  inlet  in  crosswind  in  an  upstream  irrotational  flow, 
at  low  values  of  velocity  ratio,  Uj/U,,  and/or  high  values  of  H/D,  two 
counter-rotating  vortices  trailed  from  the  rear  of  the  inlet.  These  are 
depicted  in  Figure  t(a),  which  is  taken  from  the  hydrogen  bubble  flow 
visualization  studies  in  [3].  When  H/D  was  fixed  and  the  velocity  ratio 
was  increased,  or  H/D  decreased  for  a  fixed  velocity  ratio,  an  inlet 
vortex/trailing  vortex  system  was  formed,  as  indicated  in  Figure  1(b). 

As  described  in  [3],  the  presence  of  this  trailing  vortex  is 
associated  with  the  variation  in  circulation  (i.e.,  circulation  round  the 
inlet)  along  the  axial  length  of  the  inlet.  A  schematic  of  the  location  of 
the  vortex  lines  in  the  flow  shown  in  Figure  1(b),  which  indicates  this,  is 
given  in  Figure  2.  The  circulation  around  curve  C2,  which  is  just 
forward  of  the  inlet  lip  and  encloses  the  inlet  vortex,  is  T2 •  Far  from 
the  lip,  around  curve  Cj ,  it  is  argued  that  the  circulation  is  small.  If 
this  is  so,  there  must  be  a  trailing  vortex  with  circulation  of  strength 
approximately  -1*2 »  ••  indicated  in  the  figure.  It  should  be  emphasized 
that  Figure  2  is  a  very  simplified  picture,  and  that  we  will  examine  the 
vorticity  field  in  more  detail  below. 


The  previous  investigations,  vhich  were  carried  out  in  a  water  tunnel, 
were  only  qualitative  in  nature.  In  particular,  they  provided  only  very 
rough  estimates  of  the  key  fluid  mechanic  parameters:  vortex  strength  and 
position.  These  are  required  for  the  determination  of  the  velocity 
non-uniformity  presented  to  the  engine  or  the  degree  to  which  particles  are 
ingested.  The  present  study  was  thus  undertaken  to  quantify  the  observed 
phenomena;  more  precisely,  to  provide  quantitative  measurements  of  the 
three-dimensional  vortex  flow.  As  reported  in  [4] ,  a  facility  for  making 
these  measurements  has  been  designed,  built,  and  calibrated.  Some  initial 
measurements  were  also  reported,  from  which  the  overall  structure  of  the 
flow  field  could  be  inferred,  and  this  was  found  to  agree  with  the 
conceptual  picture  developed  from  our  water  tunnel  observations. 

This  paper  presents  the  first  detailed  measurements  of  the  three- 
dimensional  flow  associated  with  the  inlet  vortex.  There  are  two  main 
points  that  are  addressed.  The  first  is  to  show  quantitatively  that  the 
circulation  of  the  inlet  and  trailing  vortices  are  (as  hypothesized  in  [3]) 
essentially  equal  and  opposite.  The  second  is  to  define  the  circulation 
and  the  location  (inside  the  inlet)  of  the  inlet  vortex  for  values  of  the 
pertinent  nondimensional  parameters  which  are  representative  of  those 
occurring  for  actual  aircraft,  since  this  is  a  situation  of  strong  practical 
interest.  In  addition  to  these  points,  further  information  about  the 
vortical  structure  of  these  flows,  obtained  from  water  tunnel  flow  visuali¬ 
zation,  will  also  be  described. 

EXPERIMENTAL  FACILITY  AND  EXPERIMENT  DESIGN 

All  of  the  quantitative  measurements  reported  were  conducted  in  the 
M.I.T.  Wright  Brothers  Wind  Tunnel.  A  layout  of  the  tunnel  is  shown  in 
Figure  3.  The  tunnel  is  closed  circuit,  with  an  elliptical  test  section 
3.05  m  wide,  2.3  m  high  and  4.6  m  long.  During  this  investigation  a  2.4  m 
wide  ground  plane  was  mounted  0.46  m  above  the  bottom  of  the  test  section. 
The  leading  edge  of  this  ground  plane  was  rounded  and  drooped  to  avoid  flow 
separation. 

The  mean  velocity  in  the  test  section  is  controlled  by  either  a 
variable  pitch  propeller,  which  can  induce  velocities  up  to  80  m/s,  or  an 
auxiliary  fan  (a  large  attic  fan),  which  is  suitable  for  speeds  up  to 
roughly  3  m/s.  The  propeller  and  fan  controls,  and  the  data  acquisition 
system,  are  located  in  a  control  room  adjacent  to  the  test  section. 

Further  information  is  given  by  Markham  [5]. 

The  flow  through  the  inlet  was  created  by  a  three-stage  centrifugal 
blower  located  outside  the  tunnel.  This  was  designed  with  backward  leaning 
blades  so  as  to  obtain  as  large  a  stable  flow  range  as  possible,  in  order 
to  have  flexibility  in  the  carrying  out  of  parametric  studies.  The  inlet 
face  average  velocity,  an  important  parameter  of  the  experiment,  was 
calculated  from  the  measured  pressure  drop  across  a  calibrated  perforated 
plate  in  the  inlet  blower/duct  system.  In  the  experiments,  the  maximum 


inlet  velocity  wee  135  n/e,  so  that  compressibility  effects  were  not 
Important .  Further  details  of  the  construction  and  calibration  of  the 
facility  for  inlet  vortex  flow  measurements  are  given  in  [4]  and  [6]. 

The  inlet  model  used  is  a  circular  cylinder  with  no  centerbody.  The 
inside  diameter  was  0.15  m  and  the  outside  diameter  was  0.2  m.  The  inlet 
lip  was  designed  to  have  no  separation  over  the  range  of  conditions  of 
interest,  using  the  procedure  developed  by  Stockman  and  Boles  [7] ,[8], 
which  yielded  an  inlet  lip  geometry  composed  of  elliptical  segments. 

The  nomenclature  used  to  describe  positions  on  and  inside  the  inlet  is 
that  the  theta  ( e)  coordinate  refers  to  angular  position  around  the  inlet. 
In  the  crosswind  configuration  used,  the  inlet  is  at  ninety  degrees  to  the 
mean  flow,  9  -  0  is  upstream  (at  the  "leading  edge”),  e  >  90  is  on  the  top 
of  the  inlet,  e  ■  180  is  at  the  rear,  and  9  -  270  is  at  the  bottom.  These 
conventions  apply  for  measurements  on  both  the  inner  and  outer  surface  of 

the  inlet.  The  axial  coordinate  is  measured  along  the  inlet  axis,  with  the 
zero  location  being  at  the  inlet  lip,  and  the  radial  coordinate  is  the 
radial  distance  from  the  axis  of  the  inlet. 

Although  there  is  no  difficulty  in  obtaining  kinematic  similarity 
between  the  parameters  in  the  model  test  and  in  the  actual  situation,  this 
is  not  so  for  dynamic  similarity.  In  particular,  the  large  size  of  an 
actual  aircraft  engine  precludes  exact  Reynolds  number  scaling  for 
practical  model  tests.  For  this  reason,  the  inlet  models  were  equipped 
with  1.5  mm  diameter  "trip  wires” ,  at  ±  67  degrees.  These  prevented 
laminar  separation  or  operation  in  the  transition  regime  between 
8ubcrltlcal  and  supercritical  operation,  where  the  flow  depended  strongly 
on  small  surface  finish  irregularities.  When  the  model  was  operated  with 
no  flow  through  the  inlet,  the  use  of  the  trip  wires  produced  static 
pressure  distributions  similar  to  those  round  a  circular  cylinder  at  the 
actual  (full-scale)  Reynolds  numbers  of  interest.  Description  of  the 
considerations  associated  with  the  use  of  the  trip  wires  can  be  found  in 
[4]  and  [6]. 

Two  models  were  machined  for  the  experiments,  one  of  plexiglas  and  the 
other  of  aluminum.  The  former  was  used  for  flow  visualization,  and  the 
latter  for  the  quantitative  measurements.  It  was  constructed  with  an  array 
of  forty  static  pressure  taps  on  its  surface,  as  well  as  with  other 
locations  for  instrumentation  and  traverse  inside  the  inlet.  Figure  4 
shows  the  aluminum  model  on  its  stand,  with  the  locations  of  the  static 
pressure  taps  along  the  inlet.  There  are  four  taps  at  each  location, 
equally  spaced  around  the  circumference.  During  testing,  as  will  be 
described  below,  the  inlet  is  rotated  around  its  axis,  so  that  detailed 
resolution  of  the  flow  can  be  achieved.  The  range  of  H/D  covered  by  the 
measurements  reported  was  from  1.0  to  2.0. 

The  measurements  of  inlet  vortex  strength  for  the  inlet  in  crosswind 
in  an  upstream  irrotational  freestream  were,  in  fact,  made  inside  the 


inlet.  There  were  two  reasons  for  doing  this.  First,  if  viscous  effects 
are  small,  the  circulation  around  the  inlet  vortex  will  be  constant,  not 
only  outside,  but  also  for  some  distance  inside  the  inlet.  More  precisely, 
if  there  is  a  region  in  which  the  vortex  core  is  surrounded  by  fluid  which 
has  the  freestream  stagnation  pressure,  and  if  the  circulation  is  measured 
in  this  region,  it  will  be  the  same  no  matter  where  one  measures  it.  It 
was  more  convenient,  therefore,  to  carry  out  the  measurements  inside  the 
inlet,  as  described  below.  Several  checks  were  made  on  whether  the  basic 
assumptions  did  in  fact  hold,  including  one  measurement  of  the  circulation 
around  the  inlet  vortex  at  a  location  just  below  the  inlet.  These  show 
that  the  measurements  within  the  inlet  were  adequate.  The  second  reason 
for  taking  measurements  inside  the  inlet,  which  is  a  quite  different  one, 
is  that  this  is  the  station  where  one  really  wishes  to  have  the  velocity 
data. 

Several  types  of  flow  measuring  probes  were  used  in  these  experiments. 
In  the  crosswind  tests,  three-hole  and  five-hole  probes  were  used  inside 
the  inlet  to  determine  the  magnitude  of  the  velocity  component  needed  to 
establish  the  circulation  around  the  inlet  vortex.  The  positions  of  these 
probes  are  shown  in  Figure  4.  The  measurements  were  made  by  rotating  the 
inlet  through  360  degrees  in  11.25  degree  increments.  This  yielded 
three-hole  probe,  Kiel  probe,  and  wall  static  pressure  data  at  64  evenly 
spaced  positions  around  the  circumference.  In  order  to  calibrate  the 
three-hole  probe  offset  angle  (the  angle  between  axial  and  the  probe 
normal),  the  inlet  was  run:  1)  facing  directly  upstream,  2)  at  an  H/D  »  6 
so  that  ground  effects  are  negligible,  and  3)  at  a  low  velocity  ratio. 

Five-hole  probes  were  used  to  determine  the  velocity  in  regions  of 
significant  pitch  angle.  A  prism  five-hole  probe  was  used  to  make  nearly 
diametric  traverses  across  the  inlet,  along  the  path  indicated  in  Figure  4. 

Measurements  of  the  trailing  vortex  flow  field  vere  made  at  points  on 
a  grid,  3.3  (inside)  diameters  downstream  of  the  inlet,  using  a  hemispheri¬ 
cal  five-hole  probe.  Motorized  horizontal  and  vertical  traversers  were 
used  to  position  the  probe  at  the  desired  points.  The  probe  traverse  set 
up  in  the  wind  tunnel  is  shown  in  Figure  5.  This  probe  could  be  rotated  in 
yaw  and  in  pitch  so  that  it  could  be  maintained  within  the  t  30  degree 


range  of  calibration 


The  circulation  measurements  that  were  carried  out  on  the  inlet  vortex 
at  the  locations  beneath  the  inlet,  were  also  done  using  a  single  slanted 
hot  wire,  rotated  to  the  three  different  orientations.  This  probe  was 
traversed  in  a  contour  to  enclose  the  inlet  vortex.  Before  carrying  out 
the  measurements,  smoke  flow  visualization  was  used  to  determine  that  the 
flow  angles  incident  on  the  probe  were  within  an  acceptable  range. 

Finally,  the  position  of  the  core  of  the  inlet  vortex  was  obtained 
from  smoke  flow  visualization  also.  A  plane  of  laser  light  was  used  to 
illuminate  the  desired  location  in  the  (plexiglas)  inlet.  Significant 
scattering  occurs  only  in  regions  of  smoke  concentration,  so  that  the 
vortex  core  could  be  plainly  seen.  Since  the  vortex  tended  to  wander 
slightly,  the  technique  used  to  establish  the  position  was  to  take  a  number 
of  pictures  (a  dozen,  say}  and  average  the  location. 

REGIMES  OF  VORTEX  FORMATION 

Before  describing  the  detailed  velocity  measurements,  it  is  useful  to 
give  some  description  of  the  regimes  in  which  the  inlet  vortex  is  present. 
Liu  et  al. ,  [3],  presented  some  initial  measurements  of  this,  but  in  Figure 
6  we  show  not  only  further  measurements  done  within  the  scope  of  the  present 
investigation,  but  data  of  several  other  studies  of  this  phenomenon.  The 
points  taken  using  the  present  facility  (with  two  different  ground  plane 
configurations)  are  plotted  as  solid.  The  dashed  line  between  the  "Vortex" 
and  "No  Vortex"  regimes  has  been  drawn,  somewhat  arbitrarily,  through  these 
points. 

Supplementary  information  about  the  particular  experimental  configura¬ 
tions  used  is  given  in  Table  I,  but  several  points  should  be  noted.  First, 
the  data  from  the  P&WA  and  the  GE  test  was  obtained  using  scale  models  of 
aircraft  engine  inlets.  Second,  the  angles  between  the  mean  flow  and  the 
inlet  axis  vary  from  55  to  135  degrees.  Third,  the  inlet  model  sizes  did 
not  vary  by  more  than  50%  from  that  used  in  the  present  study,  so  the 
Reynolds  numbers  were  roughly  comparable.  None  of  the  other  tests  made  use 
of  a  boundary  layer  trip,  however.  Fourth,  out  of  the  available  data,  only 
that  which  is  reported  to  have  low  values  of  velocity  shear  is  plotted. 
Finally,  the  GE  tests  were  done  using  an  axial  flow  fan  to  induce  the  flow 
through  the  inlet,  so  the  configuration  was  one  with  a  turbomachine  and 
with  a  centerbody. 

In  spite  of  these  differences,  it  can  be  seen  that  all  the  data  fall 
closely  into  a  pattern.  (The  last  point  would  not,  in  any  case,  be 
•xpected  to  make  a  difference.)  The  presence  or  absence  of  the  boundary 
layer  trip  does  not  appear  to  make  much  difference  as  far  as  the  vortex 
formation  boundary  is  concerned.  This  was  also  found  in  tests  with  the 
cylindrical  inlet  [4],  although  it  does  not  seem  evident  that  the  strength 
of  the  vortex  is  also  independent  of  this  feature.  Within  the  data  given. 
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TABLE  I  -  STUDIES  OF  BOUNDARIES  OF  INLET  VORTEX  FORMATION 


Source 


Conditions 


[4] 


0.15  m  cylindrical  inlet, 

90°  yaw,  thick  ground  boundary  layer 


present  paper 


same  as  above  except 
thin  ground  boundary  layer 


P+WA  Tests 
(Motycka  et  al, 
P+WA  TDM  2414) 


1/11  engine  scale  model, 
135°  crosswind, 

Ro*  -  250 


P+WA  Tests 


as  above,  Ro  >  100 


NGTE  experiments  0.15  m  and  0.3  m, 

Glenny  ARC  R+M  1114  inlets,  55°  crosswind  provided  by  fan 


GE  Tests  1/16  engine  scale  model, 

60°  crosswind  provided  by  fan, 
turbofan  inside  engine  model 


II 


1 


is  a  measure  of  inlet  shear,  where  n  is  normal  to  freestream 


for  the  various  angles,  it  also  appears  that  angular  orientation  is  not  a 
strong  influence,  although  it  should  be  clear  that  this  conclusion  cannot 
hold  all  the  way  to  angles  close  to  zero  (inlet  pointing  nearly  straight 
into  the  flow).  The  general  trend  is  that  low  values  of  H/D  and  high 
values  of  Ui/Ue  tend  to  promote  vortex  formation.  Although,  again, 
this  will  not  extrapolate  all  the  way  to  zero  free  stream  velocity,  the 
limited  data  the  authors  have  seen  suggests  that  the  velocity  ratios  at 
which  the  tendency  for  vortex  formation  again  decreases  (as  Uj./Ue  is 
increased)  are  too  high  to  be  of  great  relevance  to  the  practical  aspects 
of  this  problem. 

The  central  point  to  be  seen,  then,  is  that  the  present  set  of  experi¬ 
ments,  although  carried  out  on  a  simple  inlet  geometry  at  one  value  of 
inlet  yaw  angle,  gives  overall  flow  regimes  that  are  in  accord  with 
experiments  conducted  with  more  complex  (and  more  realistic)  geometries  and 
different  angles  of  yaw.  Hence  it  is  felt  that  the  results  that  appear 
below  will  not  only  be  useful  as  far  as  shedding  light  on  mechanisms,  but 
also  as  far  as  providing  information  on  the  magnitudes  of  vortices  to  be 
expected  in  situations  of  practical  concern. 

VELOCITY  MEASUREMENTS  IN  AN  INLET  IN  CROSSWIND  (Upstream  Irrotational  Flow) 

As  stated  in  the  introduction,  the  features  of  the  vortex  system 
associated  with  an  inlet  in  a  crosswind  have  previously  only  been  determined 
qualitatively,  and  the  purpose  of  this  study  was  to  obtain  quantitative 
data  on  this  topic.  Specific  points  to  be  addressed  are  the  correspondence 
between  the  circulation  around  the  inlet  vortex  and  that  around  the 
trailing  vortex,  and  the  dependence  of  vortex  circulation  and  position  on 
the  overall  flow  parameters,  Ui/Ue,  and  H/D. 

The  measure  of  the  vortex  strength  is  the  circulation.  This  can  be 
obtained  by  (numerically)  integrating  the  velocity  around  a  closed  contour 
enclosing  the  vortex.  Flow  visualization  showed  that  the  vortex  core  was 
located  near  the  bottom  of  the  inlet.  A  possible  path  that  could  be  used 
was  therefore  a  "D"  shaped  one,  comprising  a  segment  at  constant  radius, 
near  the  wall,  and  an  approximately  diametrical  part,  as  shown  in  Figure  4. 

The  location  for  the  measurements  was  chosen  with  several  constraints 
in  mind.  These  were:  1)  to  enclose  only  the  vorticity  of  the  inlet  vortex 
(in  particular,  to  avoid  the  streamwise  vorticity  which  results  from 
secondary  flow  in  the  boundary  layer),  2)  to  position  the  probes  for  the 
circumferential  part  of  the  traverse  close  enough  to  the  wall  so  that  pitch 
angles  are  small  (less  than  five  degrees)  and  three-hole  probes  can  be 
used,  and  3)  to  avoid  measuring  at  an  axial  location  where  circumferential 


velocities  due  to  potential  flow  effects  in  the  inlet  are  significant*. 

The  axial  location  chosen  was  two  diameters  (0.3  m)  from  the  lip,  since  at 
this  distance  the  potential  flow  effects  would  be  negligible.  A  turbulent 
boundary  layer  calculation  gives  a  (conservative)  estimate  of  the  boundary 
layer  thickness  to  be  6  mm  at  this  location  (displacement  thickness  0.7 
mm).  Potential  flow  calculations  of  the  flow  due  to  a  vortex  in  an  inlet 
using  the  estimates  of  vortex  strength  developed  in  [4]  suggested  that  at 
this  radial  position  the  pitch  angles  would  be  within  the  desired  range. 
Therefore,  for  the  circumferential  measurements  the  three-hole  and  Kiel 
probes  were  positioned  at  this  location.  The  contour  used  was  thus  174 
degrees  of  a  circular  arc,  with  a  radius  of  0.92  of  the  inlet  inner  radius. 

To  show  the  main  features  of  the  circumferential  velocity  due  to  the 
inlet  vortex,  the  non-dimensional  circumferential  velocity,  Uc/Ui, 
along  this  contour  is  shown  in  Figure  7  for  a  condition  of  H/D  -  1.5  and 
Ui/Ue  -45.  Although,  as  stated,  it  is  only  the  circumferential 
velocity  between  the  angles  of  191  and  17  degrees  that  is  necessary  for  the 
circulation  calculation,  the  velocity  field  over  the  rest  of  the  circumfer¬ 
ence  is  of  interest  in  developing  a  picture  of  the  general  flow  field  in 
the  inlet.  The  plot  in  Figure  7  represents  the  average  of  three  runs  at 
this  condition.  Experimental  uncertainty  for  this  and  other  circulation 
measurements  was  estimated  at  ±  10  percent. 

A  striking  result  is  the  existence  of  two  areas  of  high  circumferen¬ 
tial  velocity,  indicating  two  vortical  regions.  One  vortex,  centered 
roughly  about  280  degrees,  is  the  inlet  vortex.  Another  vortex,  however, 
which  has  not  been  previously  observed,  and  which  has  a  circulation 
opposite  to  that  of  the  inlet  vortex,  is  centered  at  approximately  135 
degrees. 

To  complement  these  velocity  measurements,  lampblack  flow  visualiza¬ 
tion  was  used  on  the  inside  surface  of  the  plexiglas  inlet.  The  results 
are  shown  as  Figures  8(a)  and  8(b).  As  can  be  seen  in  the  former,  the 
lampblack  streaks  originating  at  270  degrees  near  the  lip  shift  to  roughly 
300  degrees  at  two  diameters  from  the  lip  so  that  the  lower  vortex  moves  in 
a  clockwise  direction  as  it  proceeds  into  the  inlet.  This  is  due  to  the 
interaction  of  the  vortex  with  the  inlet  wall,  and  can  be  thought  of  as 
resulting  from  the  velocity  field  of  an  image  vortex.  The  lampblack  traces 

*  Near  the  inlet  lip,  the  flow  has  to  turn  into  the  inlet.  The 
streamline  curvature  brings  with  it  a  consequent  cross-stream  pressure 
gradient  and  hence  velocity  non-uniformity  (from  0  to  180  degrees).  At 
locations  away  from  the  lip  inside  the  inlet,  however,  the  velocity  is 
constant  over  the  annulus  (or  would  be,  if  there  were  no  inlet  vortex).  As 
one  moves  into  the  inlet  from  the  lip,  therefore,  there  must  be  significant 
circumferential  velocities  associated  with  the  redistribution  of  flow  over 
the  annulus.  Trying  to  measure  the  circulation  in  the  presence  of  such 
effects  essentially  means  trying  to  take  the  difference  of  two  large 
numbers,  and  it  is  desirable  to  avoid  such  a  situation. 


due  to  the  second  vortex  are  shown  in  Figure  8(b).  The  streaks  originate 
at  about  150  degrees  at  the  lip  and  move  in  a  counterclockwise  direction  to 
about  145  degrees  at  the  two-diameter  location,  so  that  this  vortex  has  a 
smaller  induced  effect  than  the  lower  one. 

A  plot  of  the  velocity  components  along  the  contour  of  integration  is 
shown  in  Figure  9.  The  figure  is  a  view  looking  into  the  inlet,  with  the 
freestream  velocity  coming  from  the  right.  The  groundplane  is  not  shown 
but  it  is  below  the  inlet.  Zt  can  be  seen  that  by  far  the  more  important 
part  of  the  circulation  integration  is  from  the  circumferential  velocity 
along  the  circular  arc.  This  is  in  fact  what  one  should  expect  from 
carrying  out  simple  potential  flow  calculations  of  a  vortex  (or  two 
vortices)  in  a  circular  cylinder,  since  the  vortices  are  close  to  the 
boundaries  and  the  region  over  which  the  velocities  ar<  large  will 
be  of  the  scale  of  the  separation  distance,  and  hence  rather  localized. 

The  validity  about  the  assumptions  concerning  probe  position  can  also 
be  examined.  The  requirement  that  the  probes  be  out  of  the  boundary  layer 
and  at  a  location  where  the  potential  flow  effects  of  the  inlet  have 
decayed  correspond  to  a  situation  in  which  the  stagnation  pressure  is 
constant  and  the  static  pressure  shows  only  a  small  local  depression  due  to 
the  presence  of  the  vortex.  Figure  10  shows  representative  pressure  dis¬ 
tributions  around  the  inside  of  the  inlet  (the  conditions  are  H/D  -  1.5, 
Uj/Ue  ■  45).  The  static  pressure  can  be  seen  to  be  as  desired;  however, 
there  are  two  regions  of  stagnation  pressure  defect.  One  region  is  near 
100  degrees,  outside  the  contour  used  for  the  circulation  measurement,  and 
is  thus  not  of  concern  for  these  types  of  measurements.  The  other  region, 
however,  is  near  300  degrees,  and  therefore  does  affect  the  circulation 
measurements. 

The  potential  flow  calculations  imply,  and  flow  visualization  bears 
out,  that  these  regions  of  lower  stagnation  pressure  are  due  to  the  local 
radial  velocities  induced  by  the  vortex.  In  other  words,  the  radial  flow 
induced  by  the  vortices  causes  the  boundary  layer  to  accumulate  on  one  side 
of  the  vortices;  displaced  in  the  clockwise  direction  for  the  lower  vortex, 
and  in  the  counterclockwise  direction  for  the  upper  vortex,  in  accordance 
with  the  respective  circulations  of  each. 

After  observing  this  type  of  behavior,  some  investigation  was  made  of 
the  situation  at  different  radial  positions.  It  was  concluded  that  the  6 
mm  location  would  be  adequate.  The  measured  circulation  of  the  vortex 
will,  of  course,  be  affected  due  to  the  influx  of  vortical  fluid  from  the 
boundary  layer  into  the  inside  of  the  contour.  One  can,  however,  make  a 
crude  estimate  of  the  influx  of  vorticity  and  this  leads  to  the  conclusion 
that  the  effect  is  small.  In  addition,  a  stronger  demonstration  can  be 
made  by  carrying  out  the  measurements  at  a  different  axial  location.  This 
was  done  at  a  station  two-thirds  of  a  diameter  upstream.  While  the 
measured  circulation  at  this  upstream  location  was  slightly  larger  (7 
percent),  the  differences  were  less  than  the  estimated  experimental 
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uncertainty.  Hence,  on  both  counts,  this  effect  is  not  significant  for  the 
conclusions  of  the  present  study. 

There  was,  however,  a  marked  difference  in  one  aspect  of  the  velocity 
measurements  at  the  two  axial  locations.  TMs  was  a  clockwise  shift  in  the 
angular  position  of  the  peak  circumferential  velocity,  of  about  twelve 
degrees  for  the  lower  vortex.  The  shift  was  less,  about  five  degrees,  and 
clockwise  for  the  upper  vortex.  Both  of  these  are  consistent  with  the 
pattern  seen  during  flow  visualization  as  well  as  with  potential  flow 
estimation  of  the  lateral  drift  of  a  vortex  near  a  wall.  The  other  aspect 
of  the  measurements  at  the  upstream  location  was  that  the  stagnation 
pressure  defects  were  less,  as  one  would  expect,  since  the  radial  transport 
has  had  less  time  to  act  and  the  boundary  layers  are  thinner  at  this 
location. 

Although  some  parametric  information  will  be  given  subsequently  about 
the  inlet  vortex  circulation,  at  this  point,  we  wish  to  pursue  only  the 
verification  of  the  hypothesis  concerning  the  inlet  vortex  and  trailing 
vortex  circulation.  We  thus  confine  the  present  comments  about  inlet 
vortex  circulation  to  stating  that  for  the  case  at  which  we  measured  the 
trailing  vortex  (as  described  in  the  next  section)  the  non-dimensional 
circulation  was  r/UeD  -  7. 

We  have  discussed  the  rationale  behind  making  the  measurements  of 
vortex  strength  inside  the  inlet.  As  a  further  check  on  the  utility  of 
this  procedure,  the  circulation  around  the  vortex  was  also  measured 
outside  of  the  inlet.  The  measurements  were  carried  out  at  H/D  *  1.5, 

Ui/Ue  *»  45  at  two  heights  above  the  ground:  0.33D  and  0.5D.  The  probe 
used  was  a  slanted  hot  wire,  operated  at  three  orientations. 

The  flow  in  the  area  of  interest  is  very  strongly  three-dimensional. 

To  ensure  that  probe  self-interference  would  not  occur,  smoke  flow 
visualization  was  first  used  to  give  an  indication  of  the  velocity 
direction.  The  circulation  contour  was  then  traversed  in  three  sections 
using  an  appropriate  probe  axis. 

The  results  of  the  traverse  showed  that  the  circulation  round  the 
vortex  was  approximately  O.G  of  that  measured  inside  the  inlet.  There  are 
several  possible  reasons  for  this,  such  as  (small)  background  shear  in  the 
tunnel  and  experimental  uncertainty  involved  in  the  measurement,  and  these 
are  currently  under  investigation.  The  central  point,  however,  is  that  to 
quite  a  good  approximation  the  basic  concept  of  measuring  the  vortex 
strength  at  the  location  of  interest  in  the  inlet  does  give  direct 
information  on  the  strength  of  the  vortex  between  the  ground  and  the  inlet. 

TRAILING  VORTEX  MEASUREMENT 

The  trailing  vortex  circulation  was  measured  at  a  condition  where  a 
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strong  inlet  vortex  was  present;  H/D  •  1.5,  and  Ui/u«  -  45.  The 
trailing  vortex  circulation  was  calculated  from  a  velocity  map  taken  at  a 
plane  downstream  of  the  inlet,  perpendicular  to  the  direction  of  the  mean 
flow.  Smoke  flow  visualization  showed  that  a  plane  approximately  three 
outer  diameters  (four  inner  diameters)  from  the  inlet  axis  was  suitable  for 
this.  Again  there  were  several  conflicting  constraints.  ‘It  was  desired  to 
have  a  region  of  "freestream"  flow  between  the  trailing  vortex  and  the 
ground  boundary  layer;  this  limited  how  far  back  the  traverse  plane  could 
be.  Against  this,  however,  we  wished  to  have  a  clearly  defined  trailing 
vortex  for  ease  of  measurement,  and  the  rolling  up  of  this  structure  took 
several  diameters.  After  considerable  trial  and  error,  the  three  diameter 
location  was  established  as  the  optimum  choice. 

The  measurements  were  taken  using  a  five-hole  hemispherical  probe. 

They  represent  a  twelve  second  average  of  the  measurements.  This  was 
necessary  since  time  resolved  measurements  in  the  trailing  vortex/wake 
region  showed  a  low  frequency  oscillation  of  approximately  1/6  Hz.  A 
repeatability  study  was  therefore  carried  out  using  averaging  times  of  six, 
twelve,  eighteen,  and  twenty-four  seconds.  From  this  it  was  determined 
that  twelve  seconds  would  be  adequate. 

The  cross  flow  velocity  vectors  obtained  are  shown  in  Figure  1 1 .  The 
velocity  vectors  are  shown  as  seen  looking  upstream  (towards  the  inlet). 

The  trace  of  the  inlet  is  also  given  in  the  figure.  The  region  in  the 
middle  of  the  figure  was  not  measured,  but  it  is  only  the  velocity  vectors 
on  a  contour  enclosing  the  vortex  that  is  needed. 

The  circulation  round  the  trailing  vortex  obtained  from  these  measure¬ 
ments  is  very  close  to  the  same  magnitude  as  that  found  for  the  inlet 
vortex,  namely  r/UeD  ■  7.  In  addition,  as  can  be  seen  from  a  comparison 
of  Figures  9  and  1 1 ,  the  sense  of  the  inlet  vortex  and  trailing  vortex 
circulations  is  opposite.  Thus,  as  hypothesized,  the  inlet  and  trailing 
vortex  are  opposite  in  sign  and  closely  equal  in  magnitude. 

INLET  CIRCULATION  FAR  FROM  THE  LIP 

A  central  point  of  the  above  hypothesis  is  that  the  circulation  round 
the  inlet  at  locations  far  from  the  inlet  lip  would  be  small.  This  is 
certainly  plausible,  since  the  velocities  on  the  outer  surface  of  the  inlet 
are  on  the  order  of  Ue  there.  However,  we  can  make  the  statement  more 
quantitative  using  data  on  the  static  pressure  distribution  round  the 
inlet.  This  is  plotted  in  Figure  12,  which  shows  curves  of  static  pressure 
coefficient,  Cp  (Cp  -  (P-Pref )/(Pt-Pref ) )  on  the  inlet  outer  surface 
versus  circumferential  position,  for  axial  locations  ranging  from  2.5 
diameters  to  0.17  diameters  from  the  inlet  lip.  The  conditions  are  the 
same  as  those  for  the  inlet  vortex  and  trailing  vortex  measurements,  namely 
H/D  «  1.5,  Ui/Ue  >  45.  Figure  12  can  be  regarded  as  an  extension  of 
the  data  given  by  Liu  et  al.  [3],  to  a  substantially  higher  velocity  ratio. 
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It  is  not  the  intent  of  this  paper  to  present  a  general  discussion  of 
the  features  of  the  static  pressure  distribution.  These  have  been  examined 
in  [3] ,  where  it  was  shown  that  the  static  pressure  distribution  at  both 
high  (as  pictured  here)  and  low  values  of  velocity  ratio  could  be  inter¬ 
preted  quite  well  in  terms  of  the  interaction  between  an  overall  flow  into 
the  inlet  and  the  inlet/trailing  vortex  system. 


This  applies  to  Figure  12  as  well,  but  that  is  not  the  main  reason  for 
including  the  figure.  We  do  this  to  show  the  contrast  between  the  distri¬ 
butions  near  the  lip  and  those  at  axial  locations  away  from  the  lip,  in 
particular  the  last  two  locations  at  2.0  and  2.5  D.  At  these  locations  the 
static  pressure  coefficient  has  a  magnitude  of  order  unity  and  is  roughly 
symmetric  about  9-0.  (This  can  perhaps  be  seen  only  crudely  due  to  the 
scale  in  the  figure;  it  does  become  apparent  when  one  plots  the  data  on  an 
expanded  scale. ) 


These  two  features  of  the  static  pressure  distribution  imply  that  the 


circulation  is  quite  small.  We  can  estimate  the  circulation  at  these  axial 
locations  using  the  measured  static  pressure  distribution  and  neglecting 
the  fact  that  there  is  a  small  component  of  flow  along  the  inlet.  (This 
will  lead  to  an  overestimate  of  the  circulation,  so  our  results  will  be  in 
fact  conservative.)  In  carrying  out  this  calculation  one  also  has  to 
estimate  where  the  flow  separates  at  the  rear  of  the  inlet,  but:  1 )  the 
static  pressure  distribution  is  such  that  one  can  make  this  estimate  with  a 
reasonable  degree  of  confidence,  and  2)  the  conclusion  about  the  magnitude 
of  the  circulation  is  unchanged  using  any  plausible  estimate  of  separation 
location.  Since  the  location  of  the  minimum  velocity  is  very  close  to  e  ■ 
0,  the  non-dimensional  (clockwise)  circulation  can  thus  be  estimated  as: 
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where  the  notation  8aep+t  9sep_  denotes  the  estimated  upper  and  lower 
separation  locations  at  the  rear  of  the  inlet.  For  reference,  at  the  2.5  D 
axial  station  these  are  at  8  -  101  degrees  and  237  degrees  respectively— 
the  static  pressure  coefficient  is  virtually  flat  between  these  locations, 
varying  by  less  than  ±  0.02. 


The  result  of  this  rough  calculation  is  that  the  circulation  round  the 
inlet  at  a  location  far  from  the  inlet  is  less  than  0.2  Ue0.  In  other 
words,  the  circulation  at  these  locations  is  almost  two  orders  of  magnitude 
less  than  that  of  the  inlet  vortex  (or  the  trailing  vortex),  confirming 
another  "piece"  of  the  conceptual  picture  that  we  had  developed  of  this 
phenomenon. 


EFFECTS  OF  VELOCITY  RATIO  AND  H/D  ON  VORTEX  STRENGTH  AND  POSITION 

To  investigate  the  influence  of  velocity  ratio,  Ui,  and  inlet  height 
to  diameter  ratio,  H/D,  on  vortex  strength  and  position,  a  limited 


parametric  study  was  done  of  the  velocity  field  in  the  inlet.  To  show  the 
influence  of  the  height  to  diameter  ratio,  tests  were  carried  out  at  H/D  - 
1.0,  1.25,  1.5,  1.75,  1.91  and  2.0,  keeping  a  constant  value  of  Uj/Ue  - 
45.  To  examine  the  effects  of  velocity  ratio,  tests  were  carried  out  at 
Oi/Ue  -  4,  17  and  45,  at  a  value  of  H/D  of  1.0.  (The  velocity  ratios 
of  4  and  17  were  chosen  since  data  had  been  previously  taken  [3]  on  the 
inlet  surface  static  pressure  distributions  at  these  values.)  The  velocity 
measurements  were  made  in  the  inlet  at  the  axial  station  two  diameters  from 
the  inlet  lip. 

Figure  13  shows  the  tangential  velocity  as  a  function  of  circumferen¬ 
tial  position  for  five  different  values  of  H/D:  1.0,  1.5,  1.75,  1.91,  2.0. 
To  help  figure  clarity,  the  curve  for  .H/D  -  1.25  is  not  plotted,  but  it  is 
virtually  the  same  as  that  for  H/D  -  1.0.  (This  curve,  as  well  as  other 
data  that  is  only  referred  to  here,  can  be  found  in  (9].)  The  figure 
clearly  indicates  that  the  circulation  of  the  inlet  vortex  decreases  in 
strength  as  the  height  to  diameter  ratio  is  increased.  The  change  in  shape 
of  the  curve  also  implies  a  qualitative  change  in  the  position  of  the 
vortex;  this  will  indeed  be  seen  below  to  be  the  case.  What  has  in  fact 
occurred,  as  can  be  seen  by  the  broadening  of  the  velocity  distribution,  is 
that  the  vortex  at  the  higher  values  of  H/D  has  moved  away  from  the  inlet 
wall. 

The  effect  of  decreasing  the  velocity  ratio  (at  constant  H/D)  is  shown 
in  Figure  14.  Decreasing  the  velocity  ratio  reduces  the  inlet  vortex 
strength  (this  trend  is  similar  to  that  observed  with  increasing  H/D)  but 
resulted  in  a  shift  of  the  position  of  maximum  circumferential  velocity  in 
the  opposite  direction  from  that  seen  with  increasing  H/D. 

The  overall  effect  of  these  parametric  changes  on  the  vortex  circula¬ 
tion  is  shown  in  Figure  1-5.  This  is  a  plot  of  the  non-dimensional 
circulation,  r/UeD,  obtained  by  numerically  integrating  the  circumferen¬ 
tial  velocity  component  and  the  relevant  velocity  component  along  the 
diametral  traverse,  for  all  of  the  points  at  which  data  was  obtained.  The 
points  for  velocity  ratio  of  45  are  shown  as  circles,  those  for  ratios  of. 
17  and  4  as  a  square  and  a  triangle  respectively. 

The  non-dimensionalization  used  in  the  figure  is  not  the  only  one  that 
could  have  been  adopted.  There  are  other  possibilities  using  the  inlet 
velocity  and/or  the  inlet  height.  These  were  explored,  and,  when  no  clear 
advantage  was  seen,  the  parameter  used  in  the  figure  was  adopted  because 
this  is  similar  to  the  type  of  parameter  that  one  would  use  to  characterize 
the  circulation  about  a  wing,  for  example. *  (There  is,  of  course,  the 
same  information  in  any  of  the  non-dimensionalizing  schemes,  because  one 


It  is  not  quite  the  same  since  the  diameter  used  would  be  the  outer 
diameter.  However,  since  much  of  the  data  on  inlet  vortices  has  been 
reported  using  the  inner  diameter,  we  have  kept  this  format.  The  ratio  of 
inner  diameter  to  outer  diameter  here  is  .75. 
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can  go  from  one  to  the  other  if  the  values  of  H/D  and  U^/Ue  are  given.) 

The  location  of  the  core  of  the  inlet  vortex  inside  the  inlet  was 
determined  using  smoke  introduced  into  the  vortex  near  the  ground.  The 
results,  based  on  the  photographs  that  were  taken,  are  presented  in  Figure 
16.  The  location  of  the  center  of  the  core  is  felt  to  be  depicted  fairly 
accurately,  however.  The  results  for  core  size  are  likely  to  be  somewhat 
of  an  overestimate  because  the  larger  of  the  smoke  particles  do  not  follow 
the  streamlines  very  accurately  and  tend  to  move  out  of  the  core.  An 
obvious  feature  is  the  large  change  in  the  position  of  the  vortex  between 
H/D  -  2.0  and  the  other  values  of  H/D  for  a  velocity  ratio  of  45.  This  is 
consistent  with  the  qualitative  change  in  the  nature  of  the  circumferential 
velocity  profile. 

STRUCTURE  OF  THE  VORTICITY  FIELD 

It  is  clear  that  the  inlet  vortex  flow  field  is  quite  complex. 

Because  of  this,  it  appears  to  be  helpful  to  try  to  infer  some  of  the 
overall  features  of  the  vortical  structure  that  is  associated  with  this 
flow  field.  In  doing  this,  we  emphasize  that  we  are  not  concerned  here 
with  a  detailed  description  of  the  structure  of  the  turbulence  or  of  the 
flow  in  the  vortex  core.  What  we  mean,  rather,  is  the  vortidty  field  that 
is  associated  with  the  mean  flow.  It  is  an  appreciation  of  this  latter 
that  is  the  first  step  to  an  understanding  of  the  general  phenomena  of 
interest.  The  description  that  follows  has  been  obtained  not  only  from  the 
velocity  measurements  that  have  been  presented,  but  also  from  some  flow 
visualization  experiments  that  were  carried  out  in  a  water  tunnel,  and 
these  will  be  discussed  where  appropriate. 

Effect  of  Velocity  Ratio  on  Vorticity  Field 

As  a  framework  for  describing  the  vorticity  field,  let  us  consider  the 
flow  regimes  encountered  with  an  inlet  operating  near  a  groundplane  as  the 
ratio  of  inlet  velocity  to  upstream  velocity  is  increased.  The  inlet  is  in 
crosswind  in  a  flow  that  is  uniform  far  upstream. 

With  no  inlet  velocity,  Uj/U0  ”  0,  the  inlet  is  essentially  a  semi¬ 
infinite  cylinder.  There  are  vortex  layers  shed  from  the  top  and  bottom  of 
the  inlet  and  from  the  end  of  the  inlet.  The  vortex  filaments  in  these 
layers  are  shown  schematically  in  Figure  17(a).  They  are  depicted  as 
straight  for  simplicity,  although  it  is  recognized  that  this  will  not  be 
the  case  in  the  actual  situation,  where  there  could  be  significant  deforma¬ 
tion  and  thus  formation  of  streamwise  vorticity. 

As  the  velocity  ratio  is  increased  to  a  low  value  (between  one  and 
ten,  say,  for  typical  aircraft  inlet  H/D),  the  situation  will  change. 

There  will  be  two  trailing  vortices.  The  vorticity  in  these  trailing 
vortices  is  produced  in  the  boundary  layers  on  the  outside  of  the  cylinder. 
Flow  visualization  shows  that,  near  the  inlet  lip,  the  fluid  in  the 


boundary  layara  on  tha  top  and  bottom  of  tha  inlat  la  sucked  into  the 
inlet.  Since  tha  fluid  particles  go  into  the  inlet,  and  since  vortex  lines 
are  never  cut,  the  vortex  lines  must  be  as  shown,  going  into  the  inlet. 

The  streamlines  are  sketched  at  tha  left.  The  upstream  motion  in  this 
regime  is  very  localized  near  the  lip,  and  is  not  shown. 

As  the  velocity  ratio  increases,  the  trailing  vortex  system  shifts 
counterclockwise,  as  indicated  in  Figure  17(c).  The  upper  vortex  also 
increases  in  diameter  and  the  lower  one  decreases.  The  circulation  around 
each  vortex  also  increases. 

As  the  velocity  ratio  is  increased  still  further,  the  lower  vortex 
shifts  abruptly  to  attach  to  the  ground,  forming  th*  inlet  vortex  (Figure 
17(d)). 

This  shift  was  examined  in  detail  in  water  tunnel  experiments.  The 
actual  circumstances  were  that  as  the  velocity  ratio  was  very  gradually 
increased,  the  two  trailing  vortex  structures  would  become  increasingly 
unsteady  and  the  lower  vortex  would  intermittently  attach  to  the  ground. 
Then,  as  the  velocity  ratio  was  increased  slightly  more,  a  strong,  stable, 
inlet  vortex  system  would  suddenly  appear.  Increasing  the  velocity  ratio 
still  more  resulted  mainly  in  a  decrease  in  the  core  size,  as  might  be 
expected  from  consideration  of  the  balance  between  convection  and  diffusion 
that  occurs  in  the  core. 

The  abruptness  of  the  transition  between  the  two  states  was  in  fact 
suggested  by  Marble  [10]  from  consideration  of  the  kinematics  of  trailing 
vortices  in  an  inviscid  fluid.  His  basic  argument  was  that  the  lower 
trailing  vortex  cannot  simply  shift  downwards  until  it  touched  the  ground 
at  some  far  downstream  location,  with  the  point  at  which  it  touched  then 
moving  gradually  upstream  as  velocity  ratio  increased.  Rather,  once  it  had 
touched  the  ground,  the  high  (infinite  for  a  line  vortex!)  induced 
velocities  would  result  in  a  very  rapid  motion  of  this  point,  initially  at 
right  angles  to  the  upstream  flow  and  then  upstream,  leading  to  the 
formation  of  the  inlet  vortex  configuration. 

As  stated  previously,  although  Figure  17  has  been  drawn  from  the  flow 
visualization  experiments,  it  is  to  be  regarded  only  as  a  schematic,  and 
there  is  one  point  which  should  be  mentioned  in  this  context.  The 
streamlines  in  Figure  17(d),  for  example,  are  shown  as  implying  particle 
motion  downstream,  but  this  is  not  true  for  all  the  flow  in  the  trailing 
vortex.  At  the  high  velocity  ratios,  the  low  stagnation  pressure  flow  in 
the  center  of  the  core  of  the  trailing  vortex  flows  upstream  and  into  the 
inlet  over  a  distance  that  is  an  appreciable  fraction  of  the  inlet 
diameter.  This  strong  upstream  flow  again  emphasizes  that  the  vortex 
filaments  must  thread  inside  the  inlet,  and  it  is  consistent  with  the 
existence  of  the  second  region  of  high  circumferential  velocity  that  was 
found  inside  the  inlet. 
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The  transition  from  an  inlet  vortex  to  the  two  trailing  vortices  when 
the  velocity  ratio  was  decreased  also  occurred  quite  abruptly.  Haves  first 
appeared  on  the  core  of  the  inlet  vortex  and  then,  as  the  velocity  ratio 
was  lowered  past  a  critical  point,  the  vortex  was  shed  downstream  and  the 
two  trailing  vortex  systems  appeared. 

Some  Comments  on  the  Topology  of  the  Vortex  Lines 

In  Figure  2  we  presented  a  very  simplified  picture  of  the  connection 
between  the  inlet  vortex  and  the  trailing  vortex.  This  was  developed  from 
arguments  about  the  solenoidality  of  the  vorticity  field.  He  now  return 
to  this  point  and  describe,  in  somewhat  more  detail,  the  overall  pattern 
of  the  vortex  lines  inside  of  and  on  the  inlet. 

To  do  this,  we  refer  to  Figure  18.  The  segments  of  the  vortex  lines 
at  A  form  the  inlet  vortex,  with  a  sense  of  rotation  as  indicated.  These 
vortex  filaments  thread  into  the  inlet  (as  at  B)  where  they  induce  the 
region  of  high  circumferential  velocity  near  the  bottom,  inside  the  inlet 
(e.g.,  near  300  degrees  in  Figure  7). 

The  vortex  filaments  near  C,  also  inside  the  inlet,  are  those  which 
induce  the  high  circumferential  velocities  near  the  rear,  inside  the  inlet 
(e.g.,  near  140  degrees  in  Figure  7).  Vortex  lines  cannot  end  in  the 
fluid  or  on  a  solid  stationary  surface  and  they  must  thus  either  go  to 
infinity  or  join  up  at  some  downstream  location.  For  purposes  of  this 
discussion,  it  does  not  really  matter  which,  and  we  have  thus  not  shown 
this  aspect  of  the  flow.  The  central  point  is  that  these  are  regarded  as 
having  essentially  equal  and  opposite  circulations. * 

The  vortex  lines  from  C  must  come  out  the  front  of  the  inlet.  (Note 
that  when  we  state  that  the  vortex  lines  "come  out"  of  the  inlet,  the 
reference  is  to  the  sense  of  the  vorticity  vector  and  not  to  the  fluid 
velocities.  The  latter  are,  in  fact,  into  the  inlet  in  the  region  that  we 
are  describing.)  This  is  shown  near  0  in  the  figure.  The  origin  of  this 
vorticity  must  be  in  the  boundary  layer  on  the  inlet,  since  this  is  the 
only  place  in  which  viscous  forces  are  acting.  The  vortex  filaments  are 
therefore  shown  lying  along  the  inlet.  The  flow  visualization  studies 
show  that  the  flow  separates  from  the  rear  of  the  inlet  [6],  and  that  it 

wIt  is  clear  that  because  of  the  no-slip  condition,  the  net  (area 
integral)  component  of  vorticity  in  the  axial  direction  inside  the  inlet 
must  be  zero.  This  is  not,  however,  what  we  mean.  He  are  referring  to  the 
vorticity  that  is  convected  into  the  inlet  (or  produced  during  the  convec¬ 
tion),  rather  than  that  produced  very  near  the  walls  by  viscous  forces 
acting  within  the  inlet.  For  example,  the  net  axial  vorticity  in  the 
boundary  layer  between  5  and  190  degrees  is  closely  equal  and  opposite  to 
that  in  the  inlet  vortex,  but  it  is  only  the  latter  that  we  are  discussing 


is  this  separated  flow  that  is  associated  with  the  trailing  vortex.  In 
addition,  the  static  pressure  distribution  shown  in  Figure  12  implies  that 
the  circulation  around  the  inlet  increases  over  a  distance  on  the  order  of 
one  to  several  outer  diameters  -  in  other  words  not  abruptly  on  the  length 
scale  of  a  diameter.  These  two  ideas  lead  one  to  draw  the  trailing  vortex 
structure  as  a  distribution  of  vorticity  rather  than  a  concentrated  vortex 
line.  This  is  shown  schematically  in  the  Figure  at  E,  for  the  vortex 
filaments  on  the  rear  surface  of  the  inlet  and  at  F  for  the  trailing 
vortex,  although  there  will  be  a  strong  tendency  for  the  latter  to  wrap 
roll  up  as  one  proceeds  in  a  downstream  direction. 

Figure  17  does  show  the  answers  to  several  questions  concerning  the 
topology  of  the  vortex  lines  in  and  on  the  inlet.  However,  it  by  no  means 
resolves  all  the  features  of  this  complicated  flow.  For  example,  we  have 
only  shown  the  trailing  vortex  filaments  that  are  associated  with  the  net 
circulation  about  the  trailing  vortex.  On  a  more  detailed  level,  there 
should  also  be  trailing  vorticity,  not  only  from  the  separation  line  at 
the  rear  of  the  inlet,  but  also  from  the  separation  line  that  is  present 
near  the  bottom  of  the  inlet.  The  trailing  vorticity  at  the  shed  into  the 
flow  at  this  latter  position  will  be  of  opposite  sign  than  that  shed  at 
the  former,  and  the  details  of  the  interaction  between  this  vorticity  and 
that  shed  from  the  upper  part  of  the  inlet  have  not  at  all  been  addressed 
here. 

SUMMARY  AND  CONCLUSIONS 

1 .  Detailed  measurements  and  flow  visualization  have  been  carried  out  for 
an  inlet  operating  in  crosswind.  The  purpose  was  to  obtain  quantitative 
information  on  the  structure  of  the  three-dimensional  flow  field 
associated  with  an  inlet  vortex. 

2.  The  studies  were  made  over  a  range  of  height  to  diameter  ratio,  H/D, 
and  velocity  ratio,  Ui/Ue,  so  that  the  flow  regimes  encountered  varied 
from  those  in  which  a  strong  inlet  vortex  and  a  trailing  vortex  were 
present,  to  those  in  which  no  inlet  vortex  existed,  but  two  counter¬ 
rotating  vorticies  trailed  from  the  inlet. 

3.  For  a  condition  where  a  strong  vortex  was  present,  measurements  were 
made  of  the  circulation  around  the  inlet  and  the  trailing  vortices.  It 
was  found  that  their  circulations  were  approximately  equal  and  opposite. 
This  had  been  hypothesized  previously  based  on  a  conceptual  picture  of  the 
inlet  vortex/trailing  vortex  system. 

4.  Also  for  this  condition,  the  velocity  measurements  and  the  flow 
visualization  both  showed  that  there  was  a  second  region  of  high  vorticity 
inside  the  inlet  (in  addition  to  the  inlet  vortex).  This  was  located  at 
roughly  "two  o'clock"  looking  into  the  inlet. 

5.  It  appears  that  the  axial  vortex  filaments  in  this  region  of  high 


vorticity  are  a  continuation  of  the  vortex  lines  that  are  associated  with 
the  trailing  vortex.  The  vortex  lines  in  this  region  thus  come  out  of  the 
inlet,  around  the  lip,  and  trail  downstream. 

6.  Flow  visualization  studies  were  also  made  of  the  transition  between 
the  inlet  vortex/trailing  vortex  configuration  and  the  flow  regime  with 
two  trailing  counter-rotating  vortices.  It  was  found  that  this  transition 
is  not  gradual  (as  a  function  of  Uj/Ue,  say)  but  occurs  abruptly  at  a 
critical  value  of  this  parameteer.  During  the  transition,  the  lower 
counter-rotating  or  trailing  vortex  moves  downward  to  become  the  inlet 
vortex. 

7.  A  limited  parametric  study  was  carried  out  to  quantitatively  define 
the  changes  in  vortex  strength  and  position  as  H/D  and  velocity  ratio 
Uj/Ug  are  varied.  Increasing  H/D  has  a  qualitatively,  but  not 
quantitatively,  similar  effect  to  decreasing  Uj/Ue. 
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Counter-Rotating 
Trailing  Vortices 


Figure  Is  Flow  Regimes  as  a  Function  of  Velocity  Ratio 


Figure  2:  Schematic  of  Vortex  Lines  in  Inlet  and  Trailing  Vortices 


Figure  3 :  Layout  of  MIT  Wright  Brothers  Wind  Tunnel 
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Pressure  Tap 
Locations  (9) 
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A“*"|  30cm  From  Inlet  Tip 


Traverse  System  for  Trailing  Vortex  Measurement 


Figure  6:  Boundaries  Between  Vortex  and  No  Vortex  Flow  Regimes 
(for  key  refer  to  Table  I) 


H/Ds I .50 
Uj  =  I35m/sec 


Figure  7:  Circumferential  Velocity  Distribution  Inside  Inlet 
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Figure  11:  Cross  Flow  Velocity  Vectors  in  Trailing  Vortex  Flow  Field 


D=  15  cm 

Axial  Location  :  2.5D  From  Inlet  Tip 


Figure  12 :  Static  Pressure  Distributions  on  Inlet  Surface 
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Figure  13:  Dependence  of  Circumferential  Velocity  Distribution  on  Height 
to  Diameter  Ratio,  H/D 
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Figure  16:  Vortex  Position  at  Measurement  Station 


TASK  IV:  INVESTIGATION  OF  DISCRETE-BLADE  AND  THREE-DIMENSIONAL  FLOWS 
IN  HIGHLY  LOADED  TURBOMACHINES 

A.  3D  FLOW  AND  DESIGN  METHODS  IN  HIGHLY-LOADED  TURBOMACHINES 

Considerable  progress  has  been  made  in  the  development  of  a  blade 
design  technique  for  a  3D  highly-loaded  cascade  subjected  to  sheared  inlet 
flow.  The  goal  is  an  a  priori  blade  design  allowing  for  the  effects  of 
secondary  flow.  To  date  the  treatment  is  limited  to  inviscid,  incompress¬ 
ible  flow  in  the  limit  of  very  thin  blades. 

The  method  used  combines  the  best  features  of  two  alternative  3D 
rotational  flow  descriptions:  the  "two-stream  function"  description,  and 
the  Clebsch  transformation  technique.  (In  the  former,  the  constants  of 
the  motion  are  displayed  best;  in  the  latter,  many  of  the  field  equations 
are  more  convenient. )  The  requirement  that  the  two  descriptions  yield 
identical  results  provides  useful  computational  information. 

In  both  descriptions,  the  "smoothing  technique"  (4.1), (4.2),  developed 
under  this  contract,  is  applied  for  cascades  of  large  blade  numbers  (small 
spacing).  The  result  is  a  neat,  closed  theory  (4.3)  which  allows  us  to 
determine  cascade  blade  shapes  required  to  produce  a  specified  result  in 
terms  of  cascade  performance  at  high  loading.  Numerical  examples  are  being 
prepared  to  illustrate  "good"  and  "bad"  cases.  A  computer  graphics  presen¬ 
tation  is  being  developed  to  display  the  flow  and  blade  results  in  a  manner 
to  make  them  more  understandable. 

The  method  developed  can  also  be  used  for  analysis  (performance) 


calculations 


B.  NUMERICAL  STUDIES  OF  SECONDARY  FLOW  IN  A  BEND  USING  SPECTRAL  METHODS 


Introduction 

Two  items  will  be  described  in  this  progress  report:  ( 1 )  the  computer 
code  £or  the  simulation  of  viscous  flow  in  a  bend  of  rectangular  cross- 
section  using  a  mixed  spectral  finite-difference  scheme;  and  (2)  the  formu¬ 
lation  of  spectral  element  technique  for  the  solution  of  the  Navier-Stok.es 
equation  in  the  bend  with  three  elements:  namely,  an  upstream  element,  a 
bend  element,  and  a  downstream  element. 

( 1 )  Mixed  spectral  finite  difference  code  for  viscous  flow  in  a  bend  of 
rectangular  cross-section 

This  code  has  been  designed  in  modular  form.  Each  module  performs  a 
particular  task  or  computes  a  desired  quantity.  The  strategy  followed  is 
that  we  first  compute  a  basic  flow  with  Vor  -  0,  Vo0  -  Vo0(r,z), 

Voz  -  0  satisfying  the  non-slip  boundary  condition  at  the  inner  wall  r  - 
r^  and  the  outer  wall  r  ■  rQ  as  well  as  the  continuity  and  momentum 
equations.  This  basic  flow  V0(r,z)  is  time-independent.  The  time 
dependent  part  and  the  three-dimensional  part  are  contained  in  v  so  that 
the  resulting  velocity  field  is  described  by  V  -  VQ(r,z)  +  v(r,e,z,t). 
Likewise,  the  pressure  field  may  be  similarly  described.  This  is  done  for 
|  convenience  and  is  a  procedure  commonly  used  by  research  workers  who 

simulate  turbulence  and  investigate  transitional  phenomena  in  fluid  flows 
(Refs.  4.4  and  4.5).  Note  that  no  approximation  is  implied  by  such  a 
description.  Because  we  have  assumed  the  flow  to  be  incompressible, 


therefore  a  staggered  finite  difference  grid  system  has  to  be  used  in  the 
direction  of  the  bend  so  that  the  continuity  condition  can  be  satisfied  in 
consistence  with  the  pressure  field.  Thus  we  would  compute  the  tangential 
velocity  v8  at  points  js  -  0,1,2,...J  (the  point  js  -  0  would 
correspond  to  the  inflow  boundary  while  the  point  js  *  J  corresponds  to 
the  outflow  boundary).  However,  for  consistency  we  would  compute  the 
velocity  components  vr,  vz,  and  the  total  pressure  field  at  the 
finite  difference  grid  points  j  lying  halfway  between  two  consecutive  j8. 

The  code  consists  of  two  essential  parts: 

(A)  The  preprocessing  stage.  Here  the  resolution  of  the  simulation,  the 
geometrical  parameters  (i.e.,  the  extent  of  the  bend,  the  inner  and  outer 
radii  measured  in  units  of  half  the  bend  height),  the  Reynolds  number,  the 
inlet  shear  profile  and  the  time  step  sizes  are  chosen.  The  main  program 
PREPR03  would  call  on: 

(i)  BASIC  to  compute  the  basic  flow. 

(ii)  PREPRO  to  preprocess  the  first  fractional  step  (the  non-linear 
convective  step)  by  performing  a  LU  decomposition  of  the  tridiagonal 
matrix  operator  for  advancing  each  of  the  velocity  components  forward 
in  time. 

(iii)  PSNPRS  to  do  the  preprocessing  necessary  for  the  implementation  of 
the  Poisson  solver  for  the  intermediate  pressure  field  pt  in  the 
pressure  correction  step.  To  achieve  this,  PSNPRS  would  in  turn  call 
on: 

(a)  PMEV  to  dlagonalise  the  pressure  Poisson  operator  in  the  e-direc¬ 
tion  (i.e.,  the  bend  direction)  with  Dirichlet  boundary  condition 


it  I  ■  0  and  tha  boundary  condition  pertaining  to  the  second 
derivative  with  respect  to  e  at  d  -  eQ. 

(b)  ZNEV  to  diagonalise  the  pressure  Poisson  operator  in  the  z- 
direction  with  Neumann  boundary  conditions  at  z  -  ±1. 

(c)  COLRTH  to  compute  the  radial  collocation  operator  resulting  from 
completion  of  (a)  and  (b)  above.  It  is  subsequently  diagonalised 
incorporating  the  Neumann  boundary  condition  at  r  ■  rj  and  r  -  r0. 

(iv)  psnvel  to  perform  the  preprocessing  necessary  for  the  Poisson  solvers 

that  implement  the  viscous  correction  step.  This  PSNVEL  would  in 

turn  call  on: 

(a)  PMEVEL  to  diagnonalise  the  velocity  Poisson  operator  in  the  9- 
direction  with  Dirichlet  boundary  condition  at  0  ■  0  (inflow 
boundary)  and  Neumann  boundary  condition  at  9  ■  0O  (outflow 
boundary). 

(b)  ZNEVEL  to  diagonalise  the  velocity  Poisson  operator  in  the  z- 
direction  with  Dirichlet  boundary  conditions  at  the  walls  z  -  ± 1 . 

(c)  COLRVR  to  compute  the  radial  collocation  operator  that  results 
from  the  completion  of  (a)  and  (b)  above  for  the  radial  velocity. 

(d)  COLRVT  to  compute  the  radial  collocation  operator  that  results 
from  the  completion  of  (a)  and  (b)  above  for  the  tangential 
velocity. 

(e)  COLRVZ  to  compute  the  radial  collocation  operator  that  results 
from  the  completion  of  (a)  and  (b)  above  for  the  axial  velocity. 

The  collocation  operators  of  (c),  (d)  and  (e)  are  each  subsequently 
diagonalised  incorporating  the  Dirichlet  boundary  condition  at  the  walls 
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r  »  and  r  -  rQ.  Note  that  (c),  (d)  and  (e)  have  to  be  done  separately 
because  the  radial  part  of  the  Poisson  operators  for  the  radial,  the 
tangential  and  the  axial  velocities  are  different  due  to  the  curvature  term 
in  the  equation  of  motion. 

The  results  from  this  preprocessing  stage  are  stored  in  data  files 
for  use  in  the  simulation  stage  to  be  described  next. 

(B)  Simulation  Stage.  The  main  program  BEND  would  call  on: 

(i)  OAT AIN  to  read  in  data  (including  the  initial  condition  on  velocity 
field)  in  files  created  during  the  preprocessing  stage  into  the  core 
memory. 

(ii)  EULER  which  uses  a  simple  Euler  time-stepping  scheme  (the  time-step 
size  for  this  step  is  much  smaller  compared  to  that  specified  for 
the  fractional  stepping)  to  establish  the  values  of  velocity  at  two 
initial  time  levels  necessary  for  starting  the  fractional  time 
stepping  of  the  Navier-Stokes  equation.  It  also  establishes  the 
initial  vorticity  field  and  the  initial  M-vector  (-V  x  v  x  V)  field. 

(ill)  PRESRE  computes  the  intermediate  pressure  field  pt  that  imposes 
the  continuity  condition. 

(iv)  VISCOS  implements  the  viscous  correction  step  that  imposes  the 
non-slip  boundary  conditions  at  the  walls. 

(v)  MOWEL  to  keep  two  time  levels  of  velocity  components  at  any  instant 
This  is  necessary  for  the  two-time-level  scheme  used  here  as  well  as 
for  checking  the  approach  to  steady  state  solution  if  it  exists. 

(vi)  VOTXCY  to  compute  the  vorticity  field  in  the  flow. 

(vii)  MOVMVT  to  keep  two  time  levels  of  M-vector  at  any  instance  as 


required  by  a  two- time-level  scheme. 

(viii)  MVECTR  to  compute  the  latest  value  of  M-vector  field. 

(ix)  DV1NTN  to  compute  the  tangential  derivative  of  the  velocity  compon¬ 
ents  which  will  be  used  by  S0URC1  to  calculate  the  RHS  vector  of  the 
operator  equation  that  implements  the  convective  step.  The  sub¬ 
routine  module  MOVDTN  would  perform  a  similar  function  to  MOWEL  and 
MOVMVT  by  keeping  two  time  levels  of  the  tangential  derivative  of 
velocity  components  at  any  one  instant. 

(x)  (a)  CONVEC  to  Implement  the  nonlinear  convective  step. 

(b)  PRESRE  to  implement  the  pressure  correction  step. 

(c)  VISCOS  to  implement  the  viscous  correction  step. 

The  PRESRE  would  call  on  PRSSLN  which  would  in  turn  invoke  BDYPSR , 
PSPH1Z ,  PSTMTN  and  SOLPSN  for  the  solution  at  the  pressure  step. 

Similarly,  VISCOS  would  call  on: 

(a)  VRDSLN  which  would  in  turn  invoke  BDYVRD,  VRPHIZ,  VRTMIN  and  SOLVRD 
to  obtain  the  radial  velocity  field. 

(b)  VTDSLN  which  would  in  turn  invoke  BDYVTN,  VTPHIZ,  VTTMTN  and  SOLVTN 
to  obtain  the  tangential  velocity  field. 

(c)  VZDSLN  which  would  in  turn  invoke  BDYVZ,  VZPHIZ,  VZTMTN  and  SOLVZ  to 
obtain  the  axial  velocity  field. 

Note  that  the  tangential  velocity  corresponds  to  the  velocity  in  the 
direction  of  the  bend  while  the  radial  and  axial  velocity  constitute  the 
flow  on  the  cross-sectional  plane  of  the  bend.  Steps  (v)  to  (x)  would  be 
repeated  for  advancing  the  flow  forward  in  time. 

At  the  moment,  the  subroutine  module  PRESRE  (for  the  pressure 
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correction  step)  has  not  been  performing  well  in  terms  of  accuracy  as  one 
would  expect.  Currently  work  is  directed  at  ascertaining  the  cause  for  the 
loss  in  accuracy  in  PRESRE  and  a  way  of  correcting  that. 

The  inflow  boundary  conditions  are  specified  while  the  outflow 
boundary  conditions  assume  either  the  first  derivative  or  the  second 
derivative  in  the  tangential  direction  vanishes. 

(2)  Formulation  of  Spectral  Element  Method 

Work  has  begun  on  the  formulation  of  spectral  element  technique  for 
the  simulation  of  three-dimensional  viscous  flow.  The  formulation 
described  here  was  carried  out  while  this  co-investigator  was  at  whittle 
Laboratory  for  a  period  of  5  weeks.  For  the  flow  in  a  bend,  three  macro¬ 
elements  are  used:  the  element  upstream  of  the  bend,  the  bend  element,  and 
the  element  downstream  of  the  bend.  Flow  variables  are  expanded  in  a 
triple  Chebyshev  series  (i.e.,  the  Chebyshev  series  is  used  as  a  high  order 
Lagrangian  interpolant) .  Unwanted  high  resolution  at  the  inflow  and 
outflow  boundaries  are  eliminated  via  the  introduction  of  multiple  elements 
in  the  flow  direction.  A  Cartesian  co-ordinate  system  is  used  in  the 
upstream  and  the  downstream  element  while  a  cylindrical  co-ordinate  system 
is  used  in  the  bend  element.  Initial  effort  is  focused  on  the  solution  of 
Poisson  and  Helmholtz  equations  (as  required  for  the  pressure  and  viscous 
steps  in  the  solution  of  Navier-Stokes  equation)  in  the  multielements 
domain.  Each  elemental  domain  is  transformed  to  a  standard  domain  in  which 
the  boundaries  of  the  domain  are  at  -1  or  +1.  Diagonalisation  of  the 
Poisson  and  the  Helmholtz  operators  is  performed  in  two  orthogonal 


directions  on  the  cross-sectional  plane.  This  reduces  the  original  three- 


dimensional  equation  into  a  one-dimensional  one  in  the  direction  of  the 
bend.  A  variational  principle  in  the  bend-direction  (or  e-direction)  is 
then  invoked  for  building  the  elemental  matrices.  These  elemental  matrices 
are  then  synthesized  into  a  global  matrix  that  satisfies  the  specified 
inflow  and  outflow  boundary  conditions  and  preserves  the  continuity  of 
function  at  the  elemental  interfaces.  It  is  hoped  that  in  the  future  we 
would  be  able  to  continue  on  this  work  so  that  an  executable  formulation 
and  analysis  could  be  used  to  simulate  viscous  internal  flow. 


C.  THREE-DIMENSIONAL  BLADE  DESIGN  FOR  LARGE  DEFLECTION 

Since  publishing  the  work  on  the  Cambridge  design  method  (reference 
4.7),  we  can  report  some  progress  both  on  the  application  of  the  method  and 
on  its  improvement.  The  method  has  been  applied  in  a  simpler,  but  the  same 
fundamental,  form  to  the  design  of  a  radial  flow  turbine  by  an  industrial 
firm.  The  results  were  very  encouraging.  Work  is  in  progress  on  the 
preparation  of  more  rigorous  methods  for  design- in  turbomachines  of  arbi¬ 
trary  geometry,  but  with  restriction  to  incompressible,  irrotational  flow 
through  the  blades. 

In  our  solution  method  (reference  4.8)  we  have  divided  the  flow  velo¬ 
city  into  a  mean  velocity,  averaged  in  the  pitchwise  or  tangential  direction 
and  a  periodic  velocity.  This  periodic  velocity  is  then  expressed  as  a 
Fourier  series  in  the  e  or  tangential  direction.  In  the  computation,  the 
Fourier  series  was  truncated  after  some  fifteen  terms.  The  errors  due  to 
this  truncation  have  been  examined  for  the  two-dimensional  cascade  and 
found  to  arise  from  a  term  in  the  Fourier  series  which  varies  as  1/n2 
where  n  is  the  harmonic.  This  term  has  been  identified  and  a  relatively 
simple  method  of  correcting  the  largest  part  of  the  truncation  error  is  now 
proposed.  When  applied  to  the  two-dimensional  cascade  for  which  the  Biot 
Savart  method  gives  the  result  of  taking  an  infinite  number  of  terms,  it  is 
found  that  using  10  terms  plus  the  truncation  error  correction  gives  a  more 
accurate  answer  than  using  60  terms  without  truncation  error  correction. 

The  method  has  been  extended  to  the  full  three-dimensional  case  and  similar 
results  have  been  obtained.  A  report  is  in  preparation. 

Some  numerical  work  has  been  done  to  examine  in  detail  the  blade 


shapes  obtained  at  and  near  the  wall.  It  has  been  found  that  the  behavior 
of  the  inclination  of  the  blade  surface  to  the  radial  direction  is 
singular  owing  to  its  reflection  in  the  wall. 

Some  progress  has  also  been  made  on  the  analytical  part  of  this  work. 
In  particular,  a  theoretical  expression  was  obtained  for  the  pressure  jump 
(or  more  appropriately  the  enthalpy  jump  in  the  case  of  irrotational 
compressible  flow)  across  a  blade  of  finite  thickness  represented  by  two 
vortex  sheets  (fig.  1):  one  pertains  to  the  suction  surface  defined  by 
a"  ■  0  -  f (r,z)  -  constant  while  the  other  pertaius  to  the  pressure 
surface  defined  by  a+  ■  0  -  f+(r,z)  -  constant.  With  the  thickness 
angle  of  the  blade  given  by  r  -  a+  -  a”  ■  f“  -  f+,  the  enthalpy 
jump  (h+  -  h“)  across  the  blade  is  derived  as 


h*  - h"  •  r-  t  +  a*>  •Mo*; 


7a+  •  V(rV  )+ 
_ 0 
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where  B  is  the  number  of  blades,  Vf  is  the  relative  velocity  as  seen  from 
a  frame  fixed  to  the  blade  row  and  rV0  -  (rV9)+  +  (rV0)“. 

The  (rV0)±  are  related  to  the  W*  and  a-  as  follows: 


.  ,  Va+  •  7 (rV  )+ 

W  - -1  - — — ^ —  7a  -  7(rV  ) 
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(see  Ref.  4.4) 

It  should  be  pointed  out  that  the  expression  given  in  Eq.  (1)  does  not 
pertain  to  a  point  but  rather  to  a  fixed  z  on  the  r-0  plane  (here  we  have 
assumed  that  a  right-handed  cylindrical  co-ordinate  system  is  used).  By 
letting  x+0,  we  obtain  the  limit  of  infinitely  thin  blades  where  the 
enthalpy  jump  (h+  -  h“)  across  the  blade  is  given  by 

(h*  -  h-)  -  |1  v.bl  .  vrv9  (3) 

where  -  1/2 (W+  +  w").  This  agrees  with  the  result  obtained  in  Ref. 

4.4.  Following  this,  we  can  thus  argue  that  the  second  term  in  the 
expression  in  Eq.  (1)  is  a  correction  term  due  to  the' fact  that  the  blade 
has  a  finite  thickness.  At  this  stage,  we  have  not  as  yet  adressed  in 
rigor  the  question  of  whether  Eq.  ( 1 )  is  valid  when  the  flow  upstream  of 
the  blade  row  is  rotational.  It  should  be  kept  in  mind  that  the  above 
results  have  been  arrived  at  based  on  the  assumptions  that  the  flow  is 
inviscid,  homentropic  and  the  rothalpy  is  uniform  far  upstream. 

From  the  above  analytical  results,  one  would  notice  the  slight  increase 
in  complexity  with  the  introduction  of  thickness  into  the  blade.  In  the 
case  of  infinitely  thin  blades,  specification  of  the  swirl  schedule  (rV0) 
in  the  blade  region  is  sufficient  for  the  determination  of  the  blade  shape. 
However,  for  the  case  of  blades  with  finite  thickness,  specification  of  only 
one  of  (rV0),  (rV0)+  or  (rV0)"  is  not  sufficient.  Some  exploration 
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of  the  possibilities  has  been  done  for  the  example  of  two-dimensional 
cascades.  Strut  profiles  have  been  computed  using  the  Biot  Savart  method 
and  various  distributions  of  (rV0)+  (note  for  the  case  of  the  strut 
(rV0)+  -  (rV0)~).  For  cambered  blades  it  may  be  better  to  specify 
rV0  and  t.  For  instance,  me  could  conceive  of  the  design  calculation 
as  follows: 

(1)  Specify  rV0[»(rV0)+  +  (rV0)”]  and  the  thickness  angle  r 

(2)  Guess  an  initial  a"  (hence  a+) 

(3)  Determine  (rV0)~  (hence  (rV0)+)  approximately  from  Eq.  (2b) 
assuming  a  value  of  W“  from  an  approximate  mean  flow  calculation 

(4  Compute  the  flow  field  with  the  knowledge  of  a+,  a*,  (rV0)+  and 

<rv0r 

(5)  Update  a-  using  the  boundary  condition  on  the  suction  surface  that 
vr  •  VcT  -  0 

(6)  Update  (rV0)+  using  Eq.  (2a)  hence  (rV0)“  or  W*1  •  Va+  ■  0 

(7)  Repeat  steps  (3)  to  (6)  until  convergence  is  obtained 

Notice  that  in  the  above  procedure,  we  have  allowed  (rV0)+  and  (rV0)~ 
to  float  in  such  a  way  that  their  sum  agrees  with  the  specified  (rV0). 

A  different  procedure  from  the  above  could  have  been  followed  by 
specifying  (rV0)“  and  t  but  allowing  (rV0)  and  (rv0)+  to 
float  such  that  their  difference  is  (rV0)~.  Alternatively,  there 
may  be  the  possibility  of  specifying  the  enthalpy  distribution  and 
thickness  angle  t.  We  believe  these  many  design  procedures  could  be 
arrived  at  using  the  above  analytical  results.  We  have  not  examined  in 
detail  the  implications  of  devising  such  a  design  procedure. 


No  effort  has  been  expended  on  developing  numerical  techniques  and 
designing  a  computer  code  to  implement  any  numerical  calculations  based 
the  above  extension  of  the  theory. 
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CONCLUSIONS 


As  we  have  stated  previously,  although  the  emphasis  in  work  done 
varies  from  project  to  project,  we  do  appear  to  be  making  continued  progress 
in  all  of  the  efforts.  One  further  aspect  deserves  mention.  This  past 
year  we  had  our  first  AFRAPT  student,  Mr.  M.  Johnson.  Our  involvement  in 
the  program  has  now  expanded  and  Mr.  D.  Fink  (associated  with  General 
Electric)  and  Ms.  J.  Koo  (AVCO),  who  are  currently  students  at  MIT,  will 
also  be  starting  in  this  program  this  summer. 


